#INTRO: LOADING PACKAGES ETC:----
rm(list=ls())
cat("\014")
library(foreign)
library(reshape2)
library(data.table)
library(stargazer)
library(MatchIt)
library(ggplot2)
library(gmodels)
library(gridExtra)
library(haven)
library(plyr)
library(tidyr)
library(reshape)
library(countrycode)
library(readr)
library(gnm)
library(nls2)
library(nlstools)
library(stats)
library(dplyr)
library(matrixStats)
library(minpack.lm)
library(lubridate)
library(AER)
library(reshape)
library(reshape2)
library(xtable)
library(tikzDevice)
library(lmtest)
library(nonnest2)

# sink("./logofcode.txt")
# sink()


#FUNCTIONS:----
FE_dG=function(FE_d){ #function creating Fixed Effects
  FE_d$feC11=ifelse(FE_d$COUNTRY==11,1,0) #1 Sweden; creating country level dummies
  FE_d$feC12=ifelse(FE_d$COUNTRY==12,1,0) #2 Norway
  FE_d$feC13=ifelse(FE_d$COUNTRY==13,1,0) #3 Denmark
  FE_d$feC21=ifelse(FE_d$COUNTRY==21,1,0) #4 Belgium
  FE_d$feC22=ifelse(FE_d$COUNTRY==22,1,0) #5 Netherlands
  FE_d$feC23=ifelse(FE_d$COUNTRY==23,1,0) #6 Luxembourg
  FE_d$feC32=ifelse(FE_d$COUNTRY==32,1,0) #7 Italy
  FE_d$feC41=ifelse(FE_d$COUNTRY==41,1,0) #8 Germany
  FE_d$feC53=ifelse(FE_d$COUNTRY==53,1,0) #9 Ireland
  FEcolumns=c(paste0("feC", 1:100))#creating fixed effects matrix
  FE_m=FE_d[, names(FE_d) %in% FEcolumns]
  FE_m=as.matrix(FE_m)
}#

FE_dFSU=function(FE_d){ #function creating Fixed Effects
  FE_d$feC11=ifelse(FE_d$COUNTRY==11,1,0) #1 Sweden; creating country level dummies
  FE_d$feC12=ifelse(FE_d$COUNTRY==12,1,0) #2 Norway
  FE_d$feC13=ifelse(FE_d$COUNTRY==13,1,0) #3 Denmark
  FE_d$feC21=ifelse(FE_d$COUNTRY==21,1,0) #4 Belgium
  FE_d$feC22=ifelse(FE_d$COUNTRY==22,1,0) #5 Netherlands
  FE_d$feC23=ifelse(FE_d$COUNTRY==23,1,0) #6 Luxembourg
  FE_d$feC32=ifelse(FE_d$COUNTRY==32,1,0) #7 Italy
  FE_d$feC53=ifelse(FE_d$COUNTRY==53,1,0) #9 Ireland
  FEcolumns=c(paste0("feC", 1:100))#creating fixed effects matrix
  FE_m=FE_d[, names(FE_d) %in% FEcolumns]
  FE_m=as.matrix(FE_m)
}

NLSc=function(b0, b, X){ #function with constant
  K=((ncol(X)-41)/20)
  XX=list()
  SbX=matrix(0,nrow = nrow(X), ncol = ncol(X))
  for (k in 1:K) {
    XX[[k]]=X[1:nrow(X),(42+(k-1)*20):(41+k*20)]
    XX[[k]]=b[k]*XX[[k]]
  }
  SbX=Reduce('+', XX)
  eSbX=exp(SbX)
  eSbXG=eSbX*X[,22:41]
  eSbXGRsum=rowSums(eSbXG, na.rm = TRUE)
  eSbXGUnit=eSbXG/eSbXGRsum
  eSbXGUnitRILE=eSbXGUnit*X[,2:21]
  eSbXGUnitRILEsum=rowSums(eSbXGUnitRILE, na.rm = TRUE)
  predict=b0+eSbXGUnitRILEsum
}#

NLSfe=function(bfe, b, X, FE){ #function with fixed effects
  K=((ncol(X)-41)/20)
  XX=list()
  SbX=matrix(0,nrow = nrow(X), ncol = ncol(X))
  for (k in 1:K) {
    XX[[k]]=X[1:nrow(X),(42+(k-1)*20):(41+k*20)]
    XX[[k]]=b[k]*XX[[k]]
  }
  SbX=Reduce('+', XX)
  eSbX=exp(SbX)
  eSbXG=eSbX*X[,22:41]
  eSbXGRsum=rowSums(eSbXG, na.rm = TRUE)
  eSbXGUnit=eSbXG/eSbXGRsum
  eSbXGUnitRILE=eSbXGUnit*X[,2:21]
  eSbXGUnitRILEsum=rowSums(eSbXGUnitRILE, na.rm = TRUE)
  FE=rowSums(t((bfe*t(FE))))
  predict=FE+eSbXGUnitRILEsum
}

NLSc_dFSU1ii=function(b0, bgs, bng, Xg, Xs){ #function with constant (M7)
  K=((ncol(Xs)-41)/20)#taken Xs since Xd same length
  XXg=list()
  XXs=list()
  SbXg=matrix(0,nrow = nrow(Xg), ncol = ncol(Xg))
  SbXs=matrix(0,nrow = nrow(Xs), ncol = ncol(Xs))
  for (k in 1:K) {
    XXg[[k]]=Xg[1:nrow(Xg),(42+(k-1)*20):(41+k*20)]
    XXg[[k]]=bgs[k]*XXg[[k]]
    XXs[[k]]=Xs[1:nrow(Xs),(42+(k-1)*20):(41+k*20)]
    XXs[[k]]=bgs[k]*XXs[[k]]
  }
  SbXg=Reduce('+', XXg)
  SbXs=Reduce('+', XXs)
  eSbXg=exp(SbXg)
  eSbXs=exp(SbXs)
  eSbXgG=eSbXg*Xg[,22:41]
  eSbXsG=eSbXs*Xs[,22:41]
  eSbXsG=bng*eSbXsG#
  eSbXgGRsum=rowSums(eSbXgG, na.rm = TRUE)
  eSbXsGRsum=rowSums(eSbXsG, na.rm = TRUE)
  eSbXgGUnit=eSbXgG/(eSbXgGRsum+eSbXsGRsum)#
  eSbXsGUnit=eSbXsG/(eSbXgGRsum+eSbXsGRsum)#
  eSbXgGUnitRILE=eSbXgGUnit*Xg[,2:21]
  eSbXsGUnitRILE=eSbXsGUnit*Xs[,2:21]
  eSbXgGUnitRILEsum=rowSums(eSbXgGUnitRILE, na.rm = TRUE)
  eSbXsGUnitRILEsum=rowSums(eSbXsGUnitRILE, na.rm = TRUE)
  predict=b0+eSbXgGUnitRILEsum+eSbXsGUnitRILEsum#
}#

NLSc_dFSU1iii=function(b0, bgs, bng, bm, Xg, Xs, Xm){ #function with constant (M7)
  K=((ncol(Xs)-41)/20)#taken Xs since Xd same length
  XXg=list()
  XXs=list()
  XXm=list()
  SbXg=matrix(0,nrow = nrow(Xg), ncol = ncol(Xg))
  SbXs=matrix(0,nrow = nrow(Xs), ncol = ncol(Xs))
  SbXm=matrix(0,nrow = nrow(Xm), ncol = ncol(Xm))
  for (k in 1:K) {
    XXg[[k]]=Xg[1:nrow(Xg),(42+(k-1)*20):(41+k*20)]
    XXg[[k]]=bgs[k]*XXg[[k]]
    XXs[[k]]=Xs[1:nrow(Xs),(42+(k-1)*20):(41+k*20)]
    XXs[[k]]=bgs[k]*XXs[[k]]
    XXm[[k]]=Xm[1:nrow(Xm),(42+(k-1)*20):(41+k*20)]
    XXm[[k]]=bm[k]*XXm[[k]]
  }
  SbXg=Reduce('+', XXg)
  SbXs=Reduce('+', XXs)
  SbXm=Reduce('+', XXm)
  eSbXg=exp(SbXg)
  eSbXs=exp(SbXs)
  eSbXm=exp(SbXm)
  eSbXgG=eSbXg*Xg[,22:41]
  eSbXsG=eSbXs*Xs[,22:41]
  eSbXmG=eSbXm*Xm[,22:41]
  eSbXsG=bng*eSbXsG#
  eSbXmG=bm*eSbXmG#
  eSbXgGRsum=rowSums(eSbXgG, na.rm = TRUE)
  eSbXsGRsum=rowSums(eSbXsG, na.rm = TRUE)
  eSbXmGRsum=rowSums(eSbXmG, na.rm = TRUE)
  eSbXgGUnit=eSbXgG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXsGUnit=eSbXsG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXmGUnit=eSbXmG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXgGUnitRILE=eSbXgGUnit*Xg[,2:21]
  eSbXsGUnitRILE=eSbXsGUnit*Xs[,2:21]
  eSbXmGUnitRILE=eSbXmGUnit*Xm[,2:21]
  eSbXgGUnitRILEsum=rowSums(eSbXgGUnitRILE, na.rm = TRUE)
  eSbXsGUnitRILEsum=rowSums(eSbXsGUnitRILE, na.rm = TRUE)
  eSbXmGUnitRILEsum=rowSums(eSbXmGUnitRILE, na.rm = TRUE)
  predict=b0+eSbXgGUnitRILEsum+eSbXsGUnitRILEsum+eSbXmGUnitRILEsum#
}#

NLSfe_dFSU1ii=function(bfe, bgs, bng, Xg, Xs, FE){ #function with fixed effects (M7)
  K=((ncol(Xs)-41)/20)#taken Xs since Xd same length
  XXg=list()
  XXs=list()
  SbXg=matrix(0,nrow = nrow(Xg), ncol = ncol(Xg))
  SbXs=matrix(0,nrow = nrow(Xs), ncol = ncol(Xs))
  for (k in 1:K) {
    XXg[[k]]=Xg[1:nrow(Xg),(42+(k-1)*20):(41+k*20)]
    XXg[[k]]=bgs[k]*XXg[[k]]
    XXs[[k]]=Xs[1:nrow(Xs),(42+(k-1)*20):(41+k*20)]
    XXs[[k]]=bgs[k]*XXs[[k]]
  }
  SbXg=Reduce('+', XXg)
  SbXs=Reduce('+', XXs)
  eSbXg=exp(SbXg)
  eSbXs=exp(SbXs)
  eSbXgG=eSbXg*Xg[,22:41]
  eSbXsG=eSbXs*Xs[,22:41]
  eSbXsG=bng*eSbXsG#
  eSbXgGRsum=rowSums(eSbXgG, na.rm = TRUE)
  eSbXsGRsum=rowSums(eSbXsG, na.rm = TRUE)
  eSbXgGUnit=eSbXgG/(eSbXgGRsum+eSbXsGRsum)#
  eSbXsGUnit=eSbXsG/(eSbXgGRsum+eSbXsGRsum)#
  eSbXgGUnitRILE=eSbXgGUnit*Xg[,2:21]
  eSbXsGUnitRILE=eSbXsGUnit*Xs[,2:21]
  eSbXgGUnitRILEsum=rowSums(eSbXgGUnitRILE, na.rm = TRUE)
  eSbXsGUnitRILEsum=rowSums(eSbXsGUnitRILE, na.rm = TRUE)
  FE=rowSums(t((bfe*t(FE))))
  predict=FE+eSbXgGUnitRILEsum+eSbXsGUnitRILEsum#
}#

NLSfe_dFSU1iii=function(bfe, bgs, bng, bm, Xg, Xs, Xm, FE){ #function with fixed effects (M7)
  K=((ncol(Xs)-41)/20)#taken Xs since Xd same length
  XXg=list()
  XXs=list()
  XXm=list()
  SbXg=matrix(0,nrow = nrow(Xg), ncol = ncol(Xg))
  SbXs=matrix(0,nrow = nrow(Xs), ncol = ncol(Xs))
  SbXm=matrix(0,nrow = nrow(Xm), ncol = ncol(Xm))
  for (k in 1:K) {
    XXg[[k]]=Xg[1:nrow(Xg),(42+(k-1)*20):(41+k*20)]
    XXg[[k]]=bgs[k]*XXg[[k]]
    XXs[[k]]=Xs[1:nrow(Xs),(42+(k-1)*20):(41+k*20)]
    XXs[[k]]=bgs[k]*XXs[[k]]
    XXm[[k]]=Xm[1:nrow(Xm),(42+(k-1)*20):(41+k*20)]
    XXm[[k]]=bm[k]*XXm[[k]]
  }
  SbXg=Reduce('+', XXg)
  SbXs=Reduce('+', XXs)
  SbXm=Reduce('+', XXm)
  eSbXg=exp(SbXg)
  eSbXs=exp(SbXs)
  eSbXm=exp(SbXm)
  eSbXgG=eSbXg*Xg[,22:41]
  eSbXsG=eSbXs*Xs[,22:41]
  eSbXmG=eSbXm*Xm[,22:41]
  eSbXsG=bng*eSbXsG#
  eSbXmG=bm*eSbXmG#
  eSbXgGRsum=rowSums(eSbXgG, na.rm = TRUE)
  eSbXsGRsum=rowSums(eSbXsG, na.rm = TRUE)
  eSbXmGRsum=rowSums(eSbXmG, na.rm = TRUE)
  eSbXgGUnit=eSbXgG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXsGUnit=eSbXsG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXmGUnit=eSbXmG/(eSbXgGRsum+eSbXsGRsum+eSbXmGRsum)#
  eSbXgGUnitRILE=eSbXgGUnit*Xg[,2:21]
  eSbXsGUnitRILE=eSbXsGUnit*Xs[,2:21]
  eSbXmGUnitRILE=eSbXmGUnit*Xm[,2:21]
  eSbXgGUnitRILEsum=rowSums(eSbXgGUnitRILE, na.rm = TRUE)
  eSbXsGUnitRILEsum=rowSums(eSbXsGUnitRILE, na.rm = TRUE)
  eSbXmGUnitRILEsum=rowSums(eSbXmGUnitRILE, na.rm = TRUE)
  FE=rowSums(t((bfe*t(FE))))
  predict=FE+eSbXgGUnitRILEsum+eSbXsGUnitRILEsum+eSbXmGUnitRILEsum#
}#

#READ DATA:----
setwd("/Users/alessioalbarello/Documents/PhD/_IIYP/P24_/_dCMP")
#main dataset with CMP and Warwick's data:
dG_=read_dta("dG_final_PSRM.dta")

dG_=dG_[dG_$COUNTRY!=31,] #w/o France
dG_=dG_[dG_$Coalition==1,] #data with only Governments Coalitions
dG_WyIEy=dG_


#summary main dataset for paper:
summary(dG_$RILEg)
sd(dG_$RILEg)#15 (15.94348)
Pmed=dG_[,names(dG_) %in% paste0("RILEp", 1:20)]
Pmed=as.matrix(Pmed)
min(Pmed,na.rm=TRUE)
max(Pmed,na.rm=TRUE)
median(Pmed,na.rm=TRUE)
sd(Pmed,na.rm=TRUE)#20 (20.43498)
Padmvp=dG_[,names(dG_) %in% paste0("ADMVPp", 1:20)]
Padmvp=as.matrix(Padmvp)
mean(Padmvp,na.rm=T)#0.07384916
FORM=as.vector(unlist(c(dG_[,names(dG_) %in% c(paste0("FORMp",1:20))])))
BigPC=as.vector(unlist(c(dG_[,names(dG_) %in% c(paste0("BigPCp",1:20))])))
table(FORM,BigPC)#form is big party in coal 90/107 cases
rm(FORM,BigPC)


##subsets:
dG_WyIEn=dG_[dG_$AfterElection==1,]#data with only Governments After Elections
dG_WnIEy=dG_[dG_$W==0,] #data with only non Warwick observations
dG_WnIEn=dG_[dG_$AfterElection==1&dG_$W==0,]#data with only Governments After Elections and non Warwick observations


##read Warwick's data on Keesing’s contemporary archives' supporting parties datasets:
#FTCS dataset (declared allies parties, including formal support parties):
dS_=read_dta("dG_FTCS.dta")
dS_=dS_[dS_$COUNTRY!=31,] #w/o France
dS_=dS_[dS_$Coalition==1,] #data with only Governments Coalitions
dS_WyIEy=dS_
dS_WyIEn=dS_[dS_$AfterElection==1,]#data with only Governments After Elections
dS_WnIEy=dS_[dS_$W==0,] #data with only non Warwick observations
dS_WnIEn=dS_[dS_$AfterElection==1&dS_$W==0,]#data with only Governments After Elections and non Warwick observations
#FTCU dataset (further including undeclared allies):
dU_=read_dta("dG_FTCU.dta")
dU_=dU_[dU_$COUNTRY!=31,] #w/o France
dU_=dU_[dU_$Coalition==1,] #data with only Governments Coalitions
dU_WyIEy=dU_
dU_WyIEn=dU_[dU_$AfterElection==1,]#data with only Governments After Elections
dU_WnIEy=dU_[dU_$W==0,] #data with only non Warwick observations
dU_WnIEn=dU_[dU_$AfterElection==1&dU_$W==0,]#data with only Governments After Elections and non Warwick observations


#Table 1: seat share: dG_WyIEy, dG_WyIEyFE, dG_WyIEn, dG_WyIEnFE:----
#create matrices to pass to function:
X1columns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("SEATShCLp", 1:20))
X1_dG_WyIEy=dG_WyIEy[, names(dG_WyIEy) %in% X1columns]
X1_dG_WyIEy[is.na(X1_dG_WyIEy)]=-99
X1_dG_WyIEy=as.matrix(X1_dG_WyIEy)
X1_dG_WyIEn=dG_WyIEn[, names(dG_WyIEn) %in% X1columns]
X1_dG_WyIEn[is.na(X1_dG_WyIEn)]=-99
X1_dG_WyIEn=as.matrix(X1_dG_WyIEn)

FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
FE_dG_WyIEn=FE_dG(dG_WyIEn)

#regressions:
NLS_X1_dG_WyIEy=nls(RILEg~NLSc(b0, b, X=X1_dG_WyIEy)
                 ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEy)
summary(NLS_X1_dG_WyIEy)
R2_X1_dG_WyIEy=1-sum((residuals(NLS_X1_dG_WyIEy))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1_dG_WyIEy
R2adj_X1_dG_WyIEy=1-(1-R2_X1_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dG_WyIEy)))
R2adj_X1_dG_WyIEy
NLS_X1_dG_WyIEyFE=nls(RILEg~NLSfe(bfe, b, X=X1_dG_WyIEy, FE=FE_dG_WyIEy),start = list(b=c(0.0005)
              ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X1_dG_WyIEyFE)
R2_X1_dG_WyIEyFE=1-sum((residuals(NLS_X1_dG_WyIEyFE))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1_dG_WyIEyFE
R2adj_X1_dG_WyIEyFE=1-(1-R2_X1_dG_WyIEyFE)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dG_WyIEyFE)))
R2adj_X1_dG_WyIEyFE
NLS_X1_dG_WyIEn=nls(RILEg~NLSc(b0, b, X=X1_dG_WyIEn)
                 ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEn)
summary(NLS_X1_dG_WyIEn)
R2_X1_dG_WyIEn=1-sum((residuals(NLS_X1_dG_WyIEn))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X1_dG_WyIEn
R2adj_X1_dG_WyIEn=1-(1-R2_X1_dG_WyIEn)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dG_WyIEn)))
R2adj_X1_dG_WyIEn
NLS_X1_dG_WyIEnFE=nls(RILEg~NLSfe(bfe, b, X=X1_dG_WyIEn, FE=FE_dG_WyIEn),start = list(b=c(0.0005)
                  ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X1_dG_WyIEnFE)
R2_X1_dG_WyIEnFE=1-sum((residuals(NLS_X1_dG_WyIEnFE))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X1_dG_WyIEnFE
R2adj_X1_dG_WyIEnFE=1-(1-R2_X1_dG_WyIEnFE)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dG_WyIEnFE)))
R2adj_X1_dG_WyIEnFE

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEy)*2), ncol = 2)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp")
mod1<- lm(RILEg~ SEATShCLp, data = x)
coefs1<- NLS_X1_dG_WyIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_dG_WyIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ SEATShCLp-1, data = x)
coefs2<- NLS_X1_dG_WyIEyFE$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X1_dG_WyIEyFE)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ SEATShCLp, data = x)
coefs3<- NLS_X1_dG_WyIEn$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X1_dG_WyIEn)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ SEATShCLp-1, data = x)
coefs4<- NLS_X1_dG_WyIEnFE$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X1_dG_WyIEnFE)))
names(se4)<- names(coefs4)
stargazer(mod1, mod2, mod3, mod4,
          coef = list(coefs1, coefs2, coefs3, coefs4),
          se = list(se1, se2, se3, se4),
          covariate.labels = c("Log Seat Share"),
          type = "latex",
          out = "Table1.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_dG_WyIEy)$coefficients[2,4]
summary(NLS_X1_dG_WyIEyFE)$coefficients[1,4]
summary(NLS_X1_dG_WyIEn)$coefficients[2,4]
summary(NLS_X1_dG_WyIEnFE)$coefficients[1,4]

2*pt(abs((NLS_X1_dG_WyIEy$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEy)))[2]),summary(NLS_X1_dG_WyIEy)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WyIEyFE$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEyFE)))[1]),summary(NLS_X1_dG_WyIEyFE)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WyIEn$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEn)))[2]),summary(NLS_X1_dG_WyIEn)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WyIEnFE$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEnFE)))[1]),summary(NLS_X1_dG_WyIEnFE)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_dG_WyIEy
R2_X1_dG_WyIEyFE
R2_X1_dG_WyIEn
R2_X1_dG_WyIEnFE

#Pseudo Adjusted R2:
R2adj_X1_dG_WyIEy
R2adj_X1_dG_WyIEyFE
R2adj_X1_dG_WyIEn
R2adj_X1_dG_WyIEnFE

#Residual SE:
sigma(NLS_X1_dG_WyIEy)
sigma(NLS_X1_dG_WyIEyFE)
sigma(NLS_X1_dG_WyIEn)
sigma(NLS_X1_dG_WyIEnFE)


#Figure 1: seat share p-values:----
#create dataset for ggplot:
x=c(-500:1500)/1000
y1=2*pt(abs((0.1111935-x)/0.1685846),105,lower.tail = FALSE)
y2=2*pt(abs((0.2612711-x)/0.1870916 ),97,lower.tail = FALSE)
y3=2*pt(abs((0.1384911-x)/0.1668561 ),63,lower.tail = FALSE)
y4=2*pt(abs((0.2170318 -x)/0.1918921),55,lower.tail = FALSE)
df01=as.data.frame(cbind(x,y1))
df01$model="(1)"
colnames(df01)[colnames(df01)=="y1"]="y"
df02=as.data.frame(cbind(x,y2))
df02$model="(2)"
colnames(df02)[colnames(df02)=="y2"]="y"
df03=as.data.frame(cbind(x,y3))
df03$model="(3)"
colnames(df03)[colnames(df03)=="y3"]="y"
df04=as.data.frame(cbind(x,y4))
df04$model="(4)"
colnames(df04)[colnames(df04)=="y4"]="y"
df05=as.data.frame(cbind(x,y1))
df05$y1=0.1
df05$model="p 0.01"
colnames(df05)[colnames(df05)=="y1"]="y"
df06=as.data.frame(cbind(x,y1))
df06$y1=0.05
df06$model="p 0.05"
colnames(df06)[colnames(df06)=="y1"]="y"
df07=as.data.frame(cbind(x,y1))
df07$y1=0.01
df07$model="p 0.1"
colnames(df07)[colnames(df07)=="y1"]="y"
df2=as.data.frame(rbind(df01,df02,df03,df04,df05,df06,df07))

#plot:
plot1=ggplot(data=df2, aes(x=x, y=y, linetype=model)) +
  geom_line() +
  scale_linetype_manual(values = c("solid","twodash","dashed","6F","dotted","dotted","dotted")) +
  theme(legend.position = "top",panel.background = element_blank()
        ,panel.grid.major = element_line(size = 0.15, linetype = 'solid',colour = "grey")
        ,panel.border = element_rect(fill = "transparent",size = 0.3),
        legend.key = element_rect(fill = "transparent"),
        text = element_text(size=9))
ggsave("Figure1.pdf",plot=plot1)


#Table A.7: seat share: dG_WnIEy, dG_WnIEyFE, dG_WnIEn, dG_WnIEnFE:----
#create matrices to pass to function:
X1columns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("SEATShCLp", 1:20))
X1_dG_WnIEy=dG_WnIEy[, names(dG_WnIEy) %in% X1columns]
X1_dG_WnIEy[is.na(X1_dG_WnIEy)]=-99
X1_dG_WnIEy=as.matrix(X1_dG_WnIEy)
X1_dG_WnIEn=dG_WnIEn[, names(dG_WnIEn) %in% X1columns]
X1_dG_WnIEn[is.na(X1_dG_WnIEn)]=-99
X1_dG_WnIEn=as.matrix(X1_dG_WnIEn)

FE_dG_WnIEy=FE_dG(dG_WnIEy)#creating FE matrices
FE_dG_WnIEn=FE_dG(dG_WnIEn)

#regressions:
NLS_X1_dG_WnIEy=nls(RILEg~NLSc(b0, b, X=X1_dG_WnIEy)
                    ,start = list(b0=1,b=c(0.0005)),data=dG_WnIEy)
summary(NLS_X1_dG_WnIEy)
R2_X1_dG_WnIEy=1-sum((residuals(NLS_X1_dG_WnIEy))^2)/sum((dG_WnIEy$RILEg-mean(dG_WnIEy$RILEg))^2)
R2_X1_dG_WnIEy
R2adj_X1_dG_WnIEy=1-(1-R2_X1_dG_WnIEy)*((length(dG_WnIEy$RILEg)-1)/(df.residual(NLS_X1_dG_WnIEy)))
R2adj_X1_dG_WnIEy
NLS_X1_dG_WnIEyFE=nls(RILEg~NLSfe(bfe, b, X=X1_dG_WnIEy, FE=FE_dG_WnIEy),start = list(b=c(0.0005)
                      ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WnIEy)
summary(NLS_X1_dG_WnIEyFE)
R2_X1_dG_WnIEyFE=1-sum((residuals(NLS_X1_dG_WnIEyFE))^2)/sum((dG_WnIEy$RILEg-mean(dG_WnIEy$RILEg))^2)
R2_X1_dG_WnIEyFE
R2adj_X1_dG_WnIEyFE=1-(1-R2_X1_dG_WnIEyFE)*((length(dG_WnIEy$RILEg)-1)/(df.residual(NLS_X1_dG_WnIEyFE)))
R2adj_X1_dG_WnIEyFE
NLS_X1_dG_WnIEn=nls(RILEg~NLSc(b0, b, X=X1_dG_WnIEn)
                    ,start = list(b0=1,b=c(0.0005)),data=dG_WnIEn)
summary(NLS_X1_dG_WnIEn)
R2_X1_dG_WnIEn=1-sum((residuals(NLS_X1_dG_WnIEn))^2)/sum((dG_WnIEn$RILEg-mean(dG_WnIEn$RILEg))^2)
R2_X1_dG_WnIEn
R2adj_X1_dG_WnIEn=1-(1-R2_X1_dG_WnIEn)*((length(dG_WnIEn$RILEg)-1)/(df.residual(NLS_X1_dG_WnIEn)))
R2adj_X1_dG_WnIEn
NLS_X1_dG_WnIEnFE=nls(RILEg~NLSfe(bfe, b, X=X1_dG_WnIEn, FE=FE_dG_WnIEn),start = list(b=c(0.0005)
                      ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WnIEn)
summary(NLS_X1_dG_WnIEnFE)
R2_X1_dG_WnIEnFE=1-sum((residuals(NLS_X1_dG_WnIEnFE))^2)/sum((dG_WnIEn$RILEg-mean(dG_WnIEn$RILEg))^2)
R2_X1_dG_WnIEnFE
R2adj_X1_dG_WnIEnFE=1-(1-R2_X1_dG_WnIEnFE)*((length(dG_WnIEn$RILEg)-1)/(df.residual(NLS_X1_dG_WnIEnFE)))
R2adj_X1_dG_WnIEnFE

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WnIEy)*2), ncol = 2)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp")
mod1<- lm(RILEg~ SEATShCLp, data = x)
coefs1<- NLS_X1_dG_WnIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_dG_WnIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ SEATShCLp-1, data = x)
coefs2<- NLS_X1_dG_WnIEyFE$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X1_dG_WnIEyFE)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ SEATShCLp, data = x)
coefs3<- NLS_X1_dG_WnIEn$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X1_dG_WnIEn)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ SEATShCLp-1, data = x)
coefs4<- NLS_X1_dG_WnIEnFE$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X1_dG_WnIEnFE)))
names(se4)<- names(coefs4)
stargazer(mod1, mod2, mod3, mod4,
          coef = list(coefs1, coefs2, coefs3, coefs4),
          se = list(se1, se2, se3, se4),
          covariate.labels = c("Log Seat Share"),
          type = "latex",
          out = "TableA.7.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_dG_WnIEy)$coefficients[2,4]
summary(NLS_X1_dG_WnIEyFE)$coefficients[1,4]
summary(NLS_X1_dG_WnIEn)$coefficients[2,4]
summary(NLS_X1_dG_WnIEnFE)$coefficients[1,4]

2*pt(abs((NLS_X1_dG_WnIEy$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WnIEy)))[2]),summary(NLS_X1_dG_WnIEy)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WnIEyFE$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_dG_WnIEyFE)))[1]),summary(NLS_X1_dG_WnIEyFE)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WnIEn$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WnIEn)))[2]),summary(NLS_X1_dG_WnIEn)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X1_dG_WnIEnFE$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_dG_WnIEnFE)))[1]),summary(NLS_X1_dG_WnIEnFE)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_dG_WnIEy
R2_X1_dG_WnIEyFE
R2_X1_dG_WnIEn
R2_X1_dG_WnIEnFE

#Pseudo Adjusted R2:
R2adj_X1_dG_WnIEy
R2adj_X1_dG_WnIEyFE
R2adj_X1_dG_WnIEn
R2adj_X1_dG_WnIEnFE

#Residual SE:
sigma(NLS_X1_dG_WnIEy)
sigma(NLS_X1_dG_WnIEyFE)
sigma(NLS_X1_dG_WnIEn)
sigma(NLS_X1_dG_WnIEnFE)


#Table 2: all variables: dG_WyIEyFEn:----
#create matrices to pass to function:
X8CEDcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("DSEATShCp", 1:20), paste0("MVPCp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20)
               , paste0("BigPCp", 1:20), paste0("MCPp", 1:20), paste0("Surp", 1:20), paste0("ICPp", 1:20))
X8CED=dG_WyIEy[, names(dG_WyIEy) %in% X8CEDcolumns]
X8CED[is.na(X8CED)]=-99
X8CED=as.matrix(X8CED)
X3_134columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
X3_134=dG_WyIEy[, names(dG_WyIEy) %in% X3_134columns]
X3_134[is.na(X3_134)]=-99
X3_134=as.matrix(X3_134)
X2_14columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("FORMp", 1:20))
X2_14=dG_WyIEy[, names(dG_WyIEy) %in% X2_14columns]
X2_14[is.na(X2_14)]=-99
X2_14=as.matrix(X2_14)
X2_13columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20))
X2_13=dG_WyIEy[, names(dG_WyIEy) %in% X2_13columns]
X2_13[is.na(X2_13)]=-99
X2_13=as.matrix(X2_13)
X1_1columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20), paste0("SEATShCLp", 1:20))
X1_1=dG_WyIEy[, names(dG_WyIEy) %in% X1_1columns]
X1_1[is.na(X1_1)]=-99
X1_1=as.matrix(X1_1)

#regressions:
NLS_X8CEDc=nls(RILEg~NLSc(b0, b, X=X8CED),start = list(b0=1,b=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5))
               ,nls.control(maxiter = 200, warnOnly = T),data=dG_WyIEy)
summary(NLS_X8CEDc)
R2_X8CEDc=1-sum((residuals(NLS_X8CEDc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X8CEDc
R2adj_X8CEDc=1-(1-R2_X8CEDc)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X8CEDc)))
R2adj_X8CEDc
NLS_X3_134c=nls(RILEg~NLSc(b0, b, X=X3_134)
                ,start = list(b0=1,b=c(0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X3_134c)
R2_X3_134c=1-sum((residuals(NLS_X3_134c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X3_134c
R2adj_X3_134c=1-(1-R2_X3_134c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X3_134c)))
R2adj_X3_134c
NLS_X2_14c=nls(RILEg~NLSc(b0, b, X=X2_14)
              ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_14c)
R2_X2_14c=1-sum((residuals(NLS_X2_14c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X2_14c
R2adj_X2_14c=1-(1-R2_X2_14c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_14c)))
R2adj_X2_14c
NLS_X2_13c=nls(RILEg~NLSc(b0, b, X=X2_13)
               ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_13c)
R2_X2_13c=1-sum((residuals(NLS_X2_13c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X2_13c
R2adj_X2_13c=1-(1-R2_X2_13c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_13c)))
R2adj_X2_13c
NLS_X1_1c=nls(RILEg~NLSc(b0, b, X=X1_1)
              ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEy)
summary(NLS_X1_1c)
R2_X1_1c=1-sum((residuals(NLS_X1_1c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1_1c
R2adj_X1_1c=1-(1-R2_X1_1c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X1_1c)))
R2adj_X1_1c

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEy)*10), ncol = 10)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp", "DSEATShCLp", "MVPCp", "ADMVPCp", "FORMp", "BigPCp", "MCPp", "Surp", "ICPp")
mod8<- lm(RILEg~ SEATShCLp + DSEATShCLp + MVPCp + ADMVPCp + FORMp + BigPCp + MCPp + Surp + ICPp, data = x)
coefs8<- NLS_X8CEDc$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X8CEDc)))
names(se8)<- names(coefs8)
mod4<- lm(RILEg~ SEATShCLp + ADMVPCp + FORMp, data = x)
coefs4<- NLS_X3_134c$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X3_134c)))
names(se4)<- names(coefs4)
mod3<- lm(RILEg~ SEATShCLp + FORMp, data = x)
coefs3<- NLS_X2_14c$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X2_14c)))
names(se3)<- names(coefs3)
mod2<- lm(RILEg~ SEATShCLp + ADMVPCp, data = x)
coefs2<- NLS_X2_13c$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X2_13c)))
names(se2)<- names(coefs2)
mod1<- lm(RILEg~ SEATShCLp, data = x)
coefs1<- NLS_X1_1c$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_1c)))
names(se1)<- names(coefs1)
stargazer(mod1, mod2, mod3, mod4, mod8,
          coef = list(coefs1, coefs2, coefs3, coefs4, coefs8),
          se = list(se1, se2, se3, se4, se8),
          covariate.labels = c("Log Seat Share", "Delta Seat Share", "Median Legislative Party", "Abs Distance to MLP"
                               , "Formateur", "Biggest Coalition Party", "Median Coalition Party", "Surplus Party in the Coalition"
                               , "Inconsistent Coalition Party "),
          type = "latex",
          out = "Table2.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_1c)$coefficients[2,4]
summary(NLS_X2_13c)$coefficients[2,4]
summary(NLS_X2_14c)$coefficients[2,4]
summary(NLS_X3_134c)$coefficients[2,4]

2*pt(abs((NLS_X1_1c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_1c)))[2]),summary(NLS_X1_1c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_13c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_13c)))[2]),summary(NLS_X2_13c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_14c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_14c)))[2]),summary(NLS_X2_14c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X3_134c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X3_134c)))[2]),summary(NLS_X3_134c)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_1c
R2_X2_13c
R2_X2_14c
R2_X3_134c
R2_X8CEDc

#Pseudo Adjusted R2:
R2adj_X1_1c
R2adj_X2_13c
R2adj_X2_14c
R2adj_X3_134c
R2adj_X8CEDc

#Residual SE:
sigma(NLS_X1_1c)
sigma(NLS_X2_13c)
sigma(NLS_X2_14c)
sigma(NLS_X3_134c)
sigma(NLS_X8CEDc)


#Table A.8: all variables: dG_WyIEyFEy:----
#create matrices to pass to function:
X8CEDcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("DSEATShCp", 1:20), paste0("MVPCp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20)
               , paste0("BigPCp", 1:20), paste0("MCPp", 1:20), paste0("Surp", 1:20), paste0("ICPp", 1:20))
X8CED=dG_WyIEy[, names(dG_WyIEy) %in% X8CEDcolumns]
X8CED[is.na(X8CED)]=-99
X8CED=as.matrix(X8CED)
X3_134columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
X3_134=dG_WyIEy[, names(dG_WyIEy) %in% X3_134columns]
X3_134[is.na(X3_134)]=-99
X3_134=as.matrix(X3_134)
X2_14columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("FORMp", 1:20))
X2_14=dG_WyIEy[, names(dG_WyIEy) %in% X2_14columns]
X2_14[is.na(X2_14)]=-99
X2_14=as.matrix(X2_14)
X2_13columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20))
X2_13=dG_WyIEy[, names(dG_WyIEy) %in% X2_13columns]
X2_13[is.na(X2_13)]=-99
X2_13=as.matrix(X2_13)
X1_1columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20), paste0("SEATShCLp", 1:20))
X1_1=dG_WyIEy[, names(dG_WyIEy) %in% X1_1columns]
X1_1[is.na(X1_1)]=-99
X1_1=as.matrix(X1_1)

FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices

#regressions:
NLS_X8CEDc=nls(RILEg~NLSfe(bfe, b, X=X8CED, FE=FE_dG_WyIEy),start = list(b=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5)
               ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X8CEDc)
R2_X8CEDc=1-sum((residuals(NLS_X8CEDc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X8CEDc
R2adj_X8CEDc=1-(1-R2_X8CEDc)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X8CEDc)))
R2adj_X8CEDc
NLS_X3_134c=nls(RILEg~NLSfe(bfe, b, X=X3_134, FE=FE_dG_WyIEy),start = list(b=c(0.0005,0.0005,0.0005)
                ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X3_134c)
R2_X3_134c=1-sum((residuals(NLS_X3_134c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X3_134c
R2adj_X3_134c=1-(1-R2_X3_134c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X3_134c)))
R2adj_X3_134c
NLS_X2_14c=nls(RILEg~NLSfe(bfe, b, X=X2_14, FE=FE_dG_WyIEy),start = list(b=c(0.0005,0.0005)
                ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_14c)
R2_X2_14c=1-sum((residuals(NLS_X2_14c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X2_14c
R2adj_X2_14c=1-(1-R2_X2_14c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_14c)))
R2adj_X2_14c
NLS_X2_13c=nls(RILEg~NLSfe(bfe, b, X=X2_13, FE=FE_dG_WyIEy),start = list(b=c(0.0005,0.0005)
                ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_13c)
R2_X2_13c=1-sum((residuals(NLS_X2_13c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X2_13c
R2adj_X2_13c=1-(1-R2_X2_13c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_13c)))
R2adj_X2_13c
NLS_X1_1c=nls(RILEg~NLSfe(bfe, b, X=X1_1, FE=FE_dG_WyIEy),start = list(b=c(0.0005)
              ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X1_1c)
R2_X1_1c=1-sum((residuals(NLS_X1_1c))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1_1c
R2adj_X1_1c=1-(1-R2_X1_1c)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X1_1c)))
R2adj_X1_1c

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEy)*19), ncol = 19)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp", "DSEATShCLp", "MVPCp", "ADMVPCp", "FORMp", "BigPCp", "MCPp", "Surp", "ICPp"
                ,"FE1", "FE2", "FE3", "FE4", "FE5", "FE6", "FE7", "FE8", "FE9")
mod8<- lm(RILEg~ SEATShCLp + DSEATShCLp + MVPCp + ADMVPCp + FORMp + BigPCp + MCPp + Surp + ICPp
          + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs8<- NLS_X8CEDc$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X8CEDc)))
names(se8)<- names(coefs8)
mod4<- lm(RILEg~ SEATShCLp + ADMVPCp + FORMp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs4<- NLS_X3_134c$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X3_134c)))
names(se4)<- names(coefs4)
mod3<- lm(RILEg~ SEATShCLp + FORMp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs3<- NLS_X2_14c$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X2_14c)))
names(se3)<- names(coefs3)
mod2<- lm(RILEg~ SEATShCLp + ADMVPCp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs2<- NLS_X2_13c$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X2_13c)))
names(se2)<- names(coefs2)
mod1<- lm(RILEg~ SEATShCLp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs1<- NLS_X1_1c$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_1c)))
names(se1)<- names(coefs1)
stargazer(mod1, mod2, mod3, mod4, mod8,
          coef = list(coefs1, coefs2, coefs3, coefs4, coefs8),
          se = list(se1, se2, se3, se4, se8),
          covariate.labels = c("Log Seat Share", "Delta Seat Share", "Median Legislative Party", "Abs Distance to MLP"
                               , "Formateur", "Biggest Coalition Party", "Median Coalition Party", "Surplus Party in the Coalition"
                               , "Inconsistent Coalition Party "
                               ,"Sweden", "Norway", "Denmark", "Belgium", "Netherlands", "Luxembourg", "Italy", "Germany", "Ireland"),
          type = "latex",
          out = "TableA.8.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_1c)$coefficients[1,4]
summary(NLS_X2_13c)$coefficients[1,4]
summary(NLS_X2_14c)$coefficients[1,4]
summary(NLS_X3_134c)$coefficients[1,4]

2*pt(abs((NLS_X1_1c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_1c)))[1]),summary(NLS_X1_1c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_13c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X2_13c)))[1]),summary(NLS_X2_13c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_14c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X2_14c)))[1]),summary(NLS_X2_14c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X3_134c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X3_134c)))[1]),summary(NLS_X3_134c)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_1c
R2_X2_13c
R2_X2_14c
R2_X3_134c
R2_X8CEDc

#Pseudo Adjusted R2:
R2adj_X1_1c
R2adj_X2_13c
R2adj_X2_14c
R2adj_X3_134c
R2adj_X8CEDc

#Residual SE:
sigma(NLS_X1_1c)
sigma(NLS_X2_13c)
sigma(NLS_X2_14c)
sigma(NLS_X3_134c)
sigma(NLS_X8CEDc)


#Table A.9: all variables: dG_WyIEnFEn:----
#create matrices to pass to function:
X8CEDcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("DSEATShCp", 1:20), paste0("MVPCp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20)
               , paste0("BigPCp", 1:20), paste0("MCPp", 1:20), paste0("Surp", 1:20), paste0("ICPp", 1:20))
X8CED=dG_WyIEn[, names(dG_WyIEn) %in% X8CEDcolumns]
X8CED[is.na(X8CED)]=-99
X8CED=as.matrix(X8CED)
X3_134columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
X3_134=dG_WyIEn[, names(dG_WyIEn) %in% X3_134columns]
X3_134[is.na(X3_134)]=-99
X3_134=as.matrix(X3_134)
X2_14columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("FORMp", 1:20))
X2_14=dG_WyIEn[, names(dG_WyIEn) %in% X2_14columns]
X2_14[is.na(X2_14)]=-99
X2_14=as.matrix(X2_14)
X2_13columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20))
X2_13=dG_WyIEn[, names(dG_WyIEn) %in% X2_13columns]
X2_13[is.na(X2_13)]=-99
X2_13=as.matrix(X2_13)
X1_1columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20), paste0("SEATShCLp", 1:20))
X1_1=dG_WyIEn[, names(dG_WyIEn) %in% X1_1columns]
X1_1[is.na(X1_1)]=-99
X1_1=as.matrix(X1_1)

#regressions:
NLS_X8CEDc=nls(RILEg~NLSc(b0, b, X=X8CED),start = list(b0=1,b=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5))
               ,nls.control(maxiter = 200, warnOnly = T),data=dG_WyIEn)
summary(NLS_X8CEDc)
R2_X8CEDc=1-sum((residuals(NLS_X8CEDc))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X8CEDc
R2adj_X8CEDc=1-(1-R2_X8CEDc)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X8CEDc)))
R2adj_X8CEDc
NLS_X3_134c=nls(RILEg~NLSc(b0, b, X=X3_134)
                ,start = list(b0=1,b=c(0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X3_134c)
R2_X3_134c=1-sum((residuals(NLS_X3_134c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X3_134c
R2adj_X3_134c=1-(1-R2_X3_134c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X3_134c)))
R2adj_X3_134c
NLS_X2_14c=nls(RILEg~NLSc(b0, b, X=X2_14)
               ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X2_14c)
R2_X2_14c=1-sum((residuals(NLS_X2_14c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X2_14c
R2adj_X2_14c=1-(1-R2_X2_14c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X2_14c)))
R2adj_X2_14c
NLS_X2_13c=nls(RILEg~NLSc(b0, b, X=X2_13)
               ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X2_13c)
R2_X2_13c=1-sum((residuals(NLS_X2_13c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X2_13c
R2adj_X2_13c=1-(1-R2_X2_13c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X2_13c)))
R2adj_X2_13c
NLS_X1_1c=nls(RILEg~NLSc(b0, b, X=X1_1)
              ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEn)
summary(NLS_X1_1c)
R2_X1_1c=1-sum((residuals(NLS_X1_1c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X1_1c
R2adj_X1_1c=1-(1-R2_X1_1c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X1_1c)))
R2adj_X1_1c

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEn)*10), ncol = 10)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp", "DSEATShCLp", "MVPCp", "ADMVPCp", "FORMp", "BigPCp", "MCPp", "Surp", "ICPp")
mod8<- lm(RILEg~ SEATShCLp + DSEATShCLp + MVPCp + ADMVPCp + FORMp + BigPCp + MCPp + Surp + ICPp, data = x)
coefs8<- NLS_X8CEDc$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X8CEDc)))
names(se8)<- names(coefs8)
mod4<- lm(RILEg~ SEATShCLp + ADMVPCp + FORMp, data = x)
coefs4<- NLS_X3_134c$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X3_134c)))
names(se4)<- names(coefs4)
mod3<- lm(RILEg~ SEATShCLp + FORMp, data = x)
coefs3<- NLS_X2_14c$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X2_14c)))
names(se3)<- names(coefs3)
mod2<- lm(RILEg~ SEATShCLp + ADMVPCp, data = x)
coefs2<- NLS_X2_13c$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X2_13c)))
names(se2)<- names(coefs2)
mod1<- lm(RILEg~ SEATShCLp, data = x)
coefs1<- NLS_X1_1c$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_1c)))
names(se1)<- names(coefs1)
stargazer(mod1, mod2, mod3, mod4, mod8,
          coef = list(coefs1, coefs2, coefs3, coefs4, coefs8),
          se = list(se1, se2, se3, se4, se8),
          covariate.labels = c("Log Seat Share", "Delta Seat Share", "Median Legislative Party", "Abs Distance to MLP"
                               , "Formateur", "Biggest Coalition Party", "Median Coalition Party", "Surplus Party in the Coalition"
                               , "Inconsistent Coalition Party "),
          type = "latex",
          out = "TableA.9.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_1c)$coefficients[2,4]
summary(NLS_X2_13c)$coefficients[2,4]
summary(NLS_X2_14c)$coefficients[2,4]
summary(NLS_X3_134c)$coefficients[2,4]

2*pt(abs((NLS_X1_1c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_1c)))[2]),summary(NLS_X1_1c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_13c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_13c)))[2]),summary(NLS_X2_13c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_14c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_14c)))[2]),summary(NLS_X2_14c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X3_134c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X3_134c)))[2]),summary(NLS_X3_134c)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_1c
R2_X2_13c
R2_X2_14c
R2_X3_134c
R2_X8CEDc

#Pseudo Adjusted R2:
R2adj_X1_1c
R2adj_X2_13c
R2adj_X2_14c
R2adj_X3_134c
R2adj_X8CEDc

#Residual SE:
sigma(NLS_X1_1c)
sigma(NLS_X2_13c)
sigma(NLS_X2_14c)
sigma(NLS_X3_134c)
sigma(NLS_X8CEDc)


#Table A.10: all variables: dG_WyIEnFEy:----
#create matrices to pass to function:
X8CEDcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("DSEATShCp", 1:20), paste0("MVPCp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20)
               , paste0("BigPCp", 1:20), paste0("MCPp", 1:20), paste0("Surp", 1:20), paste0("ICPp", 1:20))
X8CED=dG_WyIEn[, names(dG_WyIEn) %in% X8CEDcolumns]
X8CED[is.na(X8CED)]=-99
X8CED=as.matrix(X8CED)
X3_134columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
X3_134=dG_WyIEn[, names(dG_WyIEn) %in% X3_134columns]
X3_134[is.na(X3_134)]=-99
X3_134=as.matrix(X3_134)
X2_14columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("FORMp", 1:20))
X2_14=dG_WyIEn[, names(dG_WyIEn) %in% X2_14columns]
X2_14[is.na(X2_14)]=-99
X2_14=as.matrix(X2_14)
X2_13columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
               , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20))
X2_13=dG_WyIEn[, names(dG_WyIEn) %in% X2_13columns]
X2_13[is.na(X2_13)]=-99
X2_13=as.matrix(X2_13)
X1_1columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20), paste0("SEATShCLp", 1:20))
X1_1=dG_WyIEn[, names(dG_WyIEn) %in% X1_1columns]
X1_1[is.na(X1_1)]=-99
X1_1=as.matrix(X1_1)

FE_dG_WyIEn=FE_dG(dG_WyIEn)#creating FE matrices

#regressions:
NLS_X8CEDc=nls(RILEg~NLSfe(bfe, b, X=X8CED, FE=FE_dG_WyIEn),start = list(b=c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5)
               ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn
               ,nls.control(maxiter = 1000, minFactor = 0.00000001, warnOnly = T))
summary(NLS_X8CEDc)
R2_X8CEDc=1-sum((residuals(NLS_X8CEDc))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X8CEDc
R2adj_X8CEDc=1-(1-R2_X8CEDc)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X8CEDc)))
R2adj_X8CEDc
NLS_X3_134c=nls(RILEg~NLSfe(bfe, b, X=X3_134, FE=FE_dG_WyIEn),start = list(b=c(0.0005,0.0005,0.0005)
                ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X3_134c)
R2_X3_134c=1-sum((residuals(NLS_X3_134c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X3_134c
R2adj_X3_134c=1-(1-R2_X3_134c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X3_134c)))
R2adj_X3_134c
NLS_X2_14c=nls(RILEg~NLSfe(bfe, b, X=X2_14, FE=FE_dG_WyIEn),start = list(b=c(0.0005,0.0005)
               ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X2_14c)
R2_X2_14c=1-sum((residuals(NLS_X2_14c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X2_14c
R2adj_X2_14c=1-(1-R2_X2_14c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X2_14c)))
R2adj_X2_14c
NLS_X2_13c=nls(RILEg~NLSfe(bfe, b, X=X2_13, FE=FE_dG_WyIEn),start = list(b=c(0.0005,0.0005)
               ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X2_13c)
R2_X2_13c=1-sum((residuals(NLS_X2_13c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X2_13c
R2adj_X2_13c=1-(1-R2_X2_13c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X2_13c)))
R2adj_X2_13c
NLS_X1_1c=nls(RILEg~NLSfe(bfe, b, X=X1_1, FE=FE_dG_WyIEn),start = list(b=c(0.0005)
              ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEn)
summary(NLS_X1_1c)
R2_X1_1c=1-sum((residuals(NLS_X1_1c))^2)/sum((dG_WyIEn$RILEg-mean(dG_WyIEn$RILEg))^2)
R2_X1_1c
R2adj_X1_1c=1-(1-R2_X1_1c)*((length(dG_WyIEn$RILEg)-1)/(df.residual(NLS_X1_1c)))
R2adj_X1_1c

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEn)*19), ncol = 19)
x<- data.frame(x)
colnames(x)<- c("RILEg", "SEATShCLp", "DSEATShCLp", "MVPCp", "ADMVPCp", "FORMp", "BigPCp", "MCPp", "Surp", "ICPp"
                ,"FE1", "FE2", "FE3", "FE4", "FE5", "FE6", "FE7", "FE8", "FE9")
mod8<- lm(RILEg~ SEATShCLp + DSEATShCLp + MVPCp + ADMVPCp + FORMp + BigPCp + MCPp + Surp + ICPp
          + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs8<- NLS_X8CEDc$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X8CEDc)))
names(se8)<- names(coefs8)
mod4<- lm(RILEg~ SEATShCLp + ADMVPCp + FORMp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs4<- NLS_X3_134c$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X3_134c)))
names(se4)<- names(coefs4)
mod3<- lm(RILEg~ SEATShCLp + FORMp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs3<- NLS_X2_14c$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X2_14c)))
names(se3)<- names(coefs3)
mod2<- lm(RILEg~ SEATShCLp + ADMVPCp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs2<- NLS_X2_13c$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X2_13c)))
names(se2)<- names(coefs2)
mod1<- lm(RILEg~ SEATShCLp + FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 + FE9 -1, data = x)
coefs1<- NLS_X1_1c$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_1c)))
names(se1)<- names(coefs1)
stargazer(mod1, mod2, mod3, mod4, mod8,
          coef = list(coefs1, coefs2, coefs3, coefs4, coefs8),
          se = list(se1, se2, se3, se4, se8),
          covariate.labels = c("Log Seat Share", "Delta Seat Share", "Median Legislative Party", "Abs Distance to MLP"
                               , "Formateur", "Biggest Coalition Party", "Median Coalition Party", "Surplus Party in the Coalition"
                               , "Inconsistent Coalition Party "
                               ,"Sweden", "Norway", "Denmark", "Belgium", "Netherlands", "Luxembourg", "Italy", "Germany", "Ireland"),
          type = "latex",
          out = "TableA.10.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_1c)$coefficients[1,4]
summary(NLS_X2_13c)$coefficients[1,4]
summary(NLS_X2_14c)$coefficients[1,4]
summary(NLS_X3_134c)$coefficients[1,4]

2*pt(abs((NLS_X1_1c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X1_1c)))[1]),summary(NLS_X1_1c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_13c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X2_13c)))[1]),summary(NLS_X2_13c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_14c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X2_14c)))[1]),summary(NLS_X2_14c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X3_134c$m$getPars()[1]-1)/sqrt(diag(vcov(NLS_X3_134c)))[1]),summary(NLS_X3_134c)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_1c
R2_X2_13c
R2_X2_14c
R2_X3_134c
R2_X8CEDc

#Pseudo Adjusted R2:
R2adj_X1_1c
R2adj_X2_13c
R2adj_X2_14c
R2adj_X3_134c
R2adj_X8CEDc

#Residual SE:
sigma(NLS_X1_1c)
sigma(NLS_X2_13c)
sigma(NLS_X2_14c)
sigma(NLS_X3_134c)
sigma(NLS_X8CEDc)

#Tables A.1-A.3: MCP: dG_WyIEy, dG_WyIEyFE, dG_WyIEn, dG_WyIEnFE:----
#create matrices to pass to function:
MCPcolumns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("SEATShCLp", 1:20), paste0("MCPp", 1:20))
MCP_dG_WyIEy=dG_WyIEy[, names(dG_WyIEy) %in% MCPcolumns]
MCP_dG_WyIEy[is.na(MCP_dG_WyIEy)]=-99
MCP_dG_WyIEy=as.matrix(MCP_dG_WyIEy)
MCP_dG_WyIEn=dG_WyIEn[, names(dG_WyIEn) %in% MCPcolumns]
MCP_dG_WyIEn[is.na(MCP_dG_WyIEn)]=-99
MCP_dG_WyIEn=as.matrix(MCP_dG_WyIEn)

#creating independent variables:
Zuwp_WyIEy=rowSums(MCP_dG_WyIEy[,2:21]*MCP_dG_WyIEy[,22:41])/rowSums(MCP_dG_WyIEy[,22:41])#un-weighted position of parties
Zuwp_WyIEn=rowSums(MCP_dG_WyIEn[,2:21]*MCP_dG_WyIEn[,22:41])/rowSums(MCP_dG_WyIEn[,22:41])
Zmcp_WyIEy=rowSums(MCP_dG_WyIEy[,2:21]*MCP_dG_WyIEy[,62:81])#position of MCP
Zmcp_WyIEn=rowSums(MCP_dG_WyIEn[,2:21]*MCP_dG_WyIEn[,62:81])

#regressions:
LM_MCP_dG_WyIEy_uwp=lm(RILEg ~ Zuwp_WyIEy, data=dG_WyIEy)#regressions with un-weighted position of parties
summary(LM_MCP_dG_WyIEy_uwp)
LM_MCP_dG_WyIEyFE_uwp=lm(RILEg ~ Zuwp_WyIEy + factor(COUNTRY) -1, data=dG_WyIEy)
summary(LM_MCP_dG_WyIEyFE_uwp)
LM_MCP_dG_WyIEn_uwp=lm(RILEg ~ Zuwp_WyIEn, data=dG_WyIEn)
names(LM_MCP_dG_WyIEn_uwp$coefficients)=c("(Intercept)","Zuwp_WyIEy")
summary(LM_MCP_dG_WyIEn_uwp)
LM_MCP_dG_WyIEnFE_uwp=lm(RILEg ~ Zuwp_WyIEn + factor(COUNTRY) -1, data=dG_WyIEn)
names(LM_MCP_dG_WyIEnFE_uwp$coefficients)=c("Zuwp_WyIEy",
      "factor(COUNTRY)11","factor(COUNTRY)12","factor(COUNTRY)13","factor(COUNTRY)21","factor(COUNTRY)22",
      "factor(COUNTRY)23","factor(COUNTRY)32","factor(COUNTRY)41","factor(COUNTRY)53")
summary(LM_MCP_dG_WyIEnFE_uwp)

LM_MCP_dG_WyIEy_mcp=lm(RILEg ~ Zmcp_WyIEy, data=dG_WyIEy)#regressions with MCP
summary(LM_MCP_dG_WyIEy_mcp)
LM_MCP_dG_WyIEyFE_mcp=lm(RILEg ~ Zmcp_WyIEy + factor(COUNTRY) -1, data=dG_WyIEy)
summary(LM_MCP_dG_WyIEyFE_mcp)
LM_MCP_dG_WyIEn_mcp=lm(RILEg ~ Zmcp_WyIEn, data=dG_WyIEn)
names(LM_MCP_dG_WyIEn_mcp$coefficients)=c("(Intercept)","Zmcp_WyIEy")
summary(LM_MCP_dG_WyIEn_mcp)
LM_MCP_dG_WyIEnFE_mcp=lm(RILEg ~ Zmcp_WyIEn + factor(COUNTRY) -1, data=dG_WyIEn)
names(LM_MCP_dG_WyIEnFE_mcp$coefficients)=c("Zmcp_WyIEy",
      "factor(COUNTRY)11","factor(COUNTRY)12","factor(COUNTRY)13","factor(COUNTRY)21","factor(COUNTRY)22",
      "factor(COUNTRY)23","factor(COUNTRY)32","factor(COUNTRY)41","factor(COUNTRY)53")
summary(LM_MCP_dG_WyIEnFE_mcp)

#exporting regressions to latex (need to manually modify output with p-values from beta=1):
#Table A.1:
stargazer(LM_MCP_dG_WyIEy_mcp,LM_MCP_dG_WyIEy_uwp
          ,covariate.labels = c("Median Coalition Party", "Un-Weighted Party Position")
          ,type = "latex",
          out = "TableA.1.tex")

#Table A.2:
LM1=LM_MCP_dG_WyIEy_mcp
LM2=LM_MCP_dG_WyIEyFE_mcp
LM3=LM_MCP_dG_WyIEn_mcp
LM4=LM_MCP_dG_WyIEnFE_mcp
stargazer(LM1,LM2,LM3,LM4,covariate.labels = c("Median Coalition Party"),type = "latex",
          out = "TableA.2.tex")

#table p-values test with null beta=1:
2*pt(abs((summary(LM_MCP_dG_WyIEy_mcp)$coefficients[2,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEy_mcp)))[2]),summary(LM_MCP_dG_WyIEy_mcp)$df[2],lower.tail = FALSE)#beta = 1
2*pt(abs((summary(LM_MCP_dG_WyIEyFE_mcp)$coefficients[1,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEyFE_mcp)))[1]),summary(LM_MCP_dG_WyIEyFE_mcp)$df[2],lower.tail = FALSE)
2*pt(abs((summary(LM_MCP_dG_WyIEn_mcp)$coefficients[2,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEn_mcp)))[2]),summary(LM_MCP_dG_WyIEn_mcp)$df[2],lower.tail = FALSE)
2*pt(abs((summary(LM_MCP_dG_WyIEnFE_mcp)$coefficients[1,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEnFE_mcp)))[1]),summary(LM_MCP_dG_WyIEnFE_mcp)$df[2],lower.tail = FALSE)

#Table A.3:
LM1=LM_MCP_dG_WyIEy_uwp
LM2=LM_MCP_dG_WyIEyFE_uwp
LM3=LM_MCP_dG_WyIEn_uwp
LM4=LM_MCP_dG_WyIEnFE_uwp
stargazer(LM1,LM2,LM3,LM4,covariate.labels = c("Un-Weighted Party Position"),type = "latex",
          out = "TableA.3.tex")
#table p-values test with null beta=1:
2*pt(abs((summary(LM_MCP_dG_WyIEy_uwp)$coefficients[2,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEy_uwp)))[2]),summary(LM_MCP_dG_WyIEy_uwp)$df[2],lower.tail = FALSE)#beta = 1
2*pt(abs((summary(LM_MCP_dG_WyIEyFE_uwp)$coefficients[1,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEyFE_uwp)))[1]),summary(LM_MCP_dG_WyIEyFE_uwp)$df[2],lower.tail = FALSE)
2*pt(abs((summary(LM_MCP_dG_WyIEn_uwp)$coefficients[2,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEn_uwp)))[2]),summary(LM_MCP_dG_WyIEn_uwp)$df[2],lower.tail = FALSE)
2*pt(abs((summary(LM_MCP_dG_WyIEnFE_uwp)$coefficients[1,1]-1)/sqrt(diag(vcov(LM_MCP_dG_WyIEnFE_uwp)))[1]),summary(LM_MCP_dG_WyIEnFE_uwp)$df[2],lower.tail = FALSE)


#Table 3: RC min parties FTCS:----

#create data for Minority parties:
dM_WyIEy=dG_#base data from dG_

dM_WyIEy$mer=paste0(dM_WyIEy$COUNTRY,dM_WyIEy$Year)#create merging variable with unique country-years
which(duplicated(dM_WyIEy$mer))
dM_WyIEy$mer[which(duplicated(dM_WyIEy$mer))]=paste0(dM_WyIEy$mer[which(duplicated(dM_WyIEy$mer))],"_")

dt=dS_WyIEy[,c("COUNTRY","Year",paste0("GOVp", 1:20,"_G"),paste0("GOVp", 1:20,"_S"))]#base data for supporting parties
dt$mer=paste0(dt$COUNTRY,dt$Year)#create merging variable with unique country-years
which(duplicated(dt$mer))
dt$mer[which(duplicated(dt$mer))]=paste0(dt$mer[which(duplicated(dt$mer))],"_")
dt=dt[,!(names(dt) %in% c("COUNTRY","Year"))]

dM_WyIEy=merge(dM_WyIEy,dt,by=c("mer"),all.x=T,all.y=F)#merge
dM_WyIEy[,c(paste0("GOVp", 1:20,"_S"))][is.na(dM_WyIEy[,c(paste0("GOVp", 1:20,"_S"))])]=0#Na to 0 for data analysis

dM_WyIEy$GOVm1=ifelse(((dM_WyIEy$GOVp1==0)&(dM_WyIEy$GOVp1_S==0)&!is.na(dM_WyIEy$RILEp1)),1,0)#dummy being in the minority
dM_WyIEy$GOVm2=ifelse(((dM_WyIEy$GOVp2==0)&(dM_WyIEy$GOVp2_S==0)&!is.na(dM_WyIEy$RILEp2)),1,0)
dM_WyIEy$GOVm3=ifelse(((dM_WyIEy$GOVp3==0)&(dM_WyIEy$GOVp3_S==0)&!is.na(dM_WyIEy$RILEp3)),1,0)
dM_WyIEy$GOVm4=ifelse(((dM_WyIEy$GOVp4==0)&(dM_WyIEy$GOVp4_S==0)&!is.na(dM_WyIEy$RILEp4)),1,0)
dM_WyIEy$GOVm5=ifelse(((dM_WyIEy$GOVp5==0)&(dM_WyIEy$GOVp5_S==0)&!is.na(dM_WyIEy$RILEp5)),1,0)
dM_WyIEy$GOVm6=ifelse(((dM_WyIEy$GOVp6==0)&(dM_WyIEy$GOVp6_S==0)&!is.na(dM_WyIEy$RILEp6)),1,0)
dM_WyIEy$GOVm7=ifelse(((dM_WyIEy$GOVp7==0)&(dM_WyIEy$GOVp7_S==0)&!is.na(dM_WyIEy$RILEp7)),1,0)
dM_WyIEy$GOVm8=ifelse(((dM_WyIEy$GOVp8==0)&(dM_WyIEy$GOVp8_S==0)&!is.na(dM_WyIEy$RILEp8)),1,0)
dM_WyIEy$GOVm9=ifelse(((dM_WyIEy$GOVp9==0)&(dM_WyIEy$GOVp9_S==0)&!is.na(dM_WyIEy$RILEp9)),1,0)
dM_WyIEy$GOVm10=ifelse(((dM_WyIEy$GOVp10==0)&(dM_WyIEy$GOVp10_S==0)&!is.na(dM_WyIEy$RILEp10)),1,0)
dM_WyIEy$GOVm11=ifelse(((dM_WyIEy$GOVp11==0)&(dM_WyIEy$GOVp11_S==0)&!is.na(dM_WyIEy$RILEp11)),1,0)
dM_WyIEy$GOVm12=ifelse(((dM_WyIEy$GOVp12==0)&(dM_WyIEy$GOVp12_S==0)&!is.na(dM_WyIEy$RILEp12)),1,0)
dM_WyIEy$GOVm13=ifelse(((dM_WyIEy$GOVp13==0)&(dM_WyIEy$GOVp13_S==0)&!is.na(dM_WyIEy$RILEp13)),1,0)
dM_WyIEy$GOVm14=ifelse(((dM_WyIEy$GOVp14==0)&(dM_WyIEy$GOVp14_S==0)&!is.na(dM_WyIEy$RILEp14)),1,0)
dM_WyIEy$GOVm15=ifelse(((dM_WyIEy$GOVp15==0)&(dM_WyIEy$GOVp15_S==0)&!is.na(dM_WyIEy$RILEp15)),1,0)
dM_WyIEy$GOVm16=ifelse(((dM_WyIEy$GOVp16==0)&(dM_WyIEy$GOVp16_S==0)&!is.na(dM_WyIEy$RILEp16)),1,0)
dM_WyIEy$GOVm17=ifelse(((dM_WyIEy$GOVp17==0)&(dM_WyIEy$GOVp17_S==0)&!is.na(dM_WyIEy$RILEp17)),1,0)
dM_WyIEy$GOVm18=ifelse(((dM_WyIEy$GOVp18==0)&(dM_WyIEy$GOVp18_S==0)&!is.na(dM_WyIEy$RILEp18)),1,0)
dM_WyIEy$GOVm19=ifelse(((dM_WyIEy$GOVp19==0)&(dM_WyIEy$GOVp19_S==0)&!is.na(dM_WyIEy$RILEp19)),1,0)
dM_WyIEy$GOVm20=ifelse(((dM_WyIEy$GOVp20==0)&(dM_WyIEy$GOVp20_S==0)&!is.na(dM_WyIEy$RILEp20)),1,0)

dM_WyIEy$st=rowSums(dM_WyIEy[, names(dM_WyIEy) %in% c(paste0("SEATSp",1:20))])#tot seats (coded)
dM_WyIEy$SEATShPLp1=log(dM_WyIEy$SEATSp1/dM_WyIEy$st)#log seat share in parliament (coded)
dM_WyIEy$SEATShPLp2=log(dM_WyIEy$SEATSp2/dM_WyIEy$st)
dM_WyIEy$SEATShPLp3=log(dM_WyIEy$SEATSp3/dM_WyIEy$st)
dM_WyIEy$SEATShPLp4=log(dM_WyIEy$SEATSp4/dM_WyIEy$st)
dM_WyIEy$SEATShPLp5=log(dM_WyIEy$SEATSp5/dM_WyIEy$st)
dM_WyIEy$SEATShPLp6=log(dM_WyIEy$SEATSp6/dM_WyIEy$st)
dM_WyIEy$SEATShPLp7=log(dM_WyIEy$SEATSp7/dM_WyIEy$st)
dM_WyIEy$SEATShPLp8=log(dM_WyIEy$SEATSp8/dM_WyIEy$st)
dM_WyIEy$SEATShPLp9=log(dM_WyIEy$SEATSp9/dM_WyIEy$st)
dM_WyIEy$SEATShPLp10=log(dM_WyIEy$SEATSp10/dM_WyIEy$st)
dM_WyIEy$SEATShPLp11=log(dM_WyIEy$SEATSp11/dM_WyIEy$st)
dM_WyIEy$SEATShPLp12=log(dM_WyIEy$SEATSp12/dM_WyIEy$st)
dM_WyIEy$SEATShPLp13=log(dM_WyIEy$SEATSp13/dM_WyIEy$st)
dM_WyIEy$SEATShPLp14=log(dM_WyIEy$SEATSp14/dM_WyIEy$st)
dM_WyIEy$SEATShPLp15=log(dM_WyIEy$SEATSp15/dM_WyIEy$st)
dM_WyIEy$SEATShPLp16=log(dM_WyIEy$SEATSp16/dM_WyIEy$st)
dM_WyIEy$SEATShPLp17=log(dM_WyIEy$SEATSp17/dM_WyIEy$st)
dM_WyIEy$SEATShPLp18=log(dM_WyIEy$SEATSp18/dM_WyIEy$st)
dM_WyIEy$SEATShPLp19=log(dM_WyIEy$SEATSp19/dM_WyIEy$st)
dM_WyIEy$SEATShPLp20=log(dM_WyIEy$SEATSp20/dM_WyIEy$st)

dM_WyIEn=dM_WyIEy[dM_WyIEy$AfterElection==1,]#data with only Governments After Elections


##average seat share governing ang non-governing parties:
G=as.matrix(dM_WyIEn[,c(paste0("GOVp",1:20))])*as.matrix(dM_WyIEn[,c(paste0("SEATShPLp",1:20))])
G[G==0|is.na(G)]=-Inf
G=exp(G)
mean(G)#avg seat share gov parties: 0.03074877
G=as.matrix(1-dM_WyIEn[,c(paste0("GOVp",1:20))])*as.matrix(dM_WyIEn[,c(paste0("SEATShPLp",1:20))])
G[G==0|is.na(G)]=-Inf
G=exp(G)
mean(G)#avg seat share non gov parties: 0.01925123
rm(G)

#create matrices to pass to function:
X1columns_g=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVp",1:20),paste0("SEATShPLp",1:20))
X1columns_s=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVp",1:20,"_S"),paste0("SEATShPLp",1:20))
X1columns_m=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVm",1:20),paste0("SEATShPLp",1:20))

X1_dM_WyIEy_g=dM_WyIEy[,c(X1columns_g)]
X1_dM_WyIEy_g[is.na(X1_dM_WyIEy_g)]=-99
X1_dM_WyIEy_g=as.matrix(X1_dM_WyIEy_g)
X1_dM_WyIEy_s=dM_WyIEy[,c(X1columns_s)]
X1_dM_WyIEy_s[is.na(X1_dM_WyIEy_s)]=-99
X1_dM_WyIEy_s=as.matrix(X1_dM_WyIEy_s)
X1_dM_WyIEy_m=dM_WyIEy[,c(X1columns_m)]
X1_dM_WyIEy_m[is.na(X1_dM_WyIEy_m)]=-99
X1_dM_WyIEy_m=as.matrix(X1_dM_WyIEy_m)

X1_dM_WyIEn_g=dM_WyIEn[,c(X1columns_g)]
X1_dM_WyIEn_g[is.na(X1_dM_WyIEn_g)]=-99
X1_dM_WyIEn_g=as.matrix(X1_dM_WyIEn_g)
X1_dM_WyIEn_s=dM_WyIEn[,c(X1columns_s)]
X1_dM_WyIEn_s[is.na(X1_dM_WyIEn_s)]=-99
X1_dM_WyIEn_s=as.matrix(X1_dM_WyIEn_s)
X1_dM_WyIEn_m=dM_WyIEn[,c(X1columns_m)]
X1_dM_WyIEn_m[is.na(X1_dM_WyIEn_m)]=-99
X1_dM_WyIEn_m=as.matrix(X1_dM_WyIEn_m)

FE_dM_WyIEy=FE_dG(dM_WyIEy)#creating FE matrices
FE_dM_WyIEn=FE_dG(dM_WyIEn)

#regressions:
NLS_X1_dS_WyIEy=nls(RILEg~NLSc_dFSU1ii(b0, bgs, bng, Xg=X1_dM_WyIEy_g, Xs=X1_dM_WyIEy_s)
                    ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005))
                    ,data=dM_WyIEy)#
summary(NLS_X1_dS_WyIEy)
nrow(dM_WyIEy)
R2_X1_dS_WyIEy=1-sum((residuals(NLS_X1_dS_WyIEy))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dS_WyIEy
R2adj_X1_dS_WyIEy=1-(1-R2_X1_dS_WyIEy)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEy)))
R2adj_X1_dS_WyIEy
NLS_X1_dSM_WyIEy=nls(RILEg~NLSc_dFSU1iii(b0,bgs,bng,bm,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,Xm=X1_dM_WyIEy_m)
                    ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005))
                    ,data=dM_WyIEy)#
summary(NLS_X1_dSM_WyIEy)
R2_X1_dSM_WyIEy=1-sum((residuals(NLS_X1_dSM_WyIEy))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dSM_WyIEy
R2adj_X1_dSM_WyIEy=1-(1-R2_X1_dSM_WyIEy)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEy)))
R2adj_X1_dSM_WyIEy

NLS_X1_dS_WyIEyFE=nls(RILEg~NLSfe_dFSU1ii(bfe,bgs,bng,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,FE=FE_dM_WyIEy)
                      ,start = list(bgs=c(0.0005),bng=c(0.0005)
                                    ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                      ,data=dM_WyIEy)
summary(NLS_X1_dS_WyIEyFE)
R2_X1_dS_WyIEyFE=1-sum((residuals(NLS_X1_dS_WyIEyFE))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dS_WyIEyFE
R2adj_X1_dS_WyIEyFE=1-(1-R2_X1_dS_WyIEyFE)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEyFE)))
R2adj_X1_dS_WyIEyFE
NLS_X1_dSM_WyIEyFE=nls(RILEg~NLSfe_dFSU1iii(bfe,bgs,bng,bm,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,Xm=X1_dM_WyIEy_m,FE=FE_dM_WyIEy)
                      ,start = list(bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005)
                                    ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                      ,data=dM_WyIEy)
summary(NLS_X1_dSM_WyIEyFE)
R2_X1_dSM_WyIEyFE=1-sum((residuals(NLS_X1_dSM_WyIEyFE))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dSM_WyIEyFE
R2adj_X1_dSM_WyIEyFE=1-(1-R2_X1_dSM_WyIEyFE)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEyFE)))
R2adj_X1_dSM_WyIEyFE

NLS_X1_dS_WyIEn=nls(RILEg~NLSc_dFSU1ii(b0, bgs, bng, Xg=X1_dM_WyIEn_g, Xs=X1_dM_WyIEn_s)
                    ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005))
                    ,data=dM_WyIEn)#
summary(NLS_X1_dS_WyIEn)
R2_X1_dS_WyIEn=1-sum((residuals(NLS_X1_dS_WyIEn))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dS_WyIEn
R2adj_X1_dS_WyIEn=1-(1-R2_X1_dS_WyIEn)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEn)))
R2adj_X1_dS_WyIEn
NLS_X1_dSM_WyIEn=nls(RILEg~NLSc_dFSU1iii(b0,bgs,bng,bm,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,Xm=X1_dM_WyIEn_m)
                     ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005))
                     ,data=dM_WyIEn)#
summary(NLS_X1_dSM_WyIEn)
R2_X1_dSM_WyIEn=1-sum((residuals(NLS_X1_dSM_WyIEn))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dSM_WyIEn
R2adj_X1_dSM_WyIEn=1-(1-R2_X1_dSM_WyIEn)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEn)))
R2adj_X1_dSM_WyIEn

NLS_X1_dS_WyIEnFE=nls(RILEg~NLSfe_dFSU1ii(bfe,bgs,bng,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,FE=FE_dM_WyIEn)
                      ,start = list(bgs=c(0.0005),bng=c(0.0005)
                                    ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                      ,data=dM_WyIEn)
summary(NLS_X1_dS_WyIEnFE)
R2_X1_dS_WyIEnFE=1-sum((residuals(NLS_X1_dS_WyIEnFE))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dS_WyIEnFE
R2adj_X1_dS_WyIEnFE=1-(1-R2_X1_dS_WyIEnFE)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEnFE)))
R2adj_X1_dS_WyIEnFE
NLS_X1_dSM_WyIEnFE=nls(RILEg~NLSfe_dFSU1iii(bfe,bgs,bng,bm,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,Xm=X1_dM_WyIEn_m,FE=FE_dM_WyIEn)
                       ,start = list(bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005)
                                     ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                       ,data=dM_WyIEn)
summary(NLS_X1_dSM_WyIEnFE)
R2_X1_dSM_WyIEnFE=1-sum((residuals(NLS_X1_dSM_WyIEnFE))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dSM_WyIEnFE
R2adj_X1_dSM_WyIEnFE=1-(1-R2_X1_dSM_WyIEnFE)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEnFE)))
R2adj_X1_dSM_WyIEnFE

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEy)*13), ncol = 13)
x<- data.frame(x)
colnames(x)<- c("RILEg", "bgs", "bng", "bm"
                ,"FE1", "FE2", "FE3", "FE4", "FE5", "FE6", "FE7", "FE8", "FE9")
mod1<- lm(RILEg~ bgs+bng, data = x)
coefs1<- NLS_X1_dS_WyIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_dS_WyIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ bgs+bng+bm, data = x)
coefs2<- NLS_X1_dSM_WyIEy$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X1_dSM_WyIEy)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ bgs+bng+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs3<- NLS_X1_dS_WyIEyFE$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X1_dS_WyIEyFE)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ bgs+bng+bm+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs4<- NLS_X1_dSM_WyIEyFE$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X1_dSM_WyIEyFE)))
names(se4)<- names(coefs4)
mod5<- lm(RILEg~ bgs+bng, data = x)
coefs5<- NLS_X1_dS_WyIEn$m$getPars()
names(coefs5)<- names(mod5$coefficients)
se5<- sqrt(diag(vcov(NLS_X1_dS_WyIEn)))
names(se5)<- names(coefs5)
mod6<- lm(RILEg~ bgs+bng+bm, data = x)
coefs6<- NLS_X1_dSM_WyIEn$m$getPars()
names(coefs6)<- names(mod6$coefficients)
se6<- sqrt(diag(vcov(NLS_X1_dSM_WyIEn)))
names(se6)<- names(coefs6)
mod7<- lm(RILEg~ bgs+bng+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs7<- NLS_X1_dS_WyIEnFE$m$getPars()
names(coefs7)<- names(mod7$coefficients)
se7<- sqrt(diag(vcov(NLS_X1_dS_WyIEnFE)))
names(se7)<- names(coefs7)
mod8<- lm(RILEg~ bgs+bng+bm+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs8<- NLS_X1_dSM_WyIEnFE$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X1_dSM_WyIEnFE)))
names(se8)<- names(coefs8)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,
          coef = list(coefs1,coefs2,coefs3,coefs4,coefs5,coefs6,coefs7,coefs8),
          se = list(se1,se2,se3,se4,se5,se6,se7,se8),
          covariate.labels = c("Log Seat Share", "Supporting Parties", "Opposition Parties"),
          type = "latex",
          out = "Table3.tex")

#Pseudo R2:
R2_X1_dS_WyIEy
R2_X1_dSM_WyIEy
R2_X1_dS_WyIEyFE
R2_X1_dSM_WyIEyFE
R2_X1_dS_WyIEn
R2_X1_dSM_WyIEn
R2_X1_dS_WyIEnFE
R2_X1_dSM_WyIEnFE

#Pseudo Adjusted R2:
R2adj_X1_dS_WyIEy
R2adj_X1_dSM_WyIEy
R2adj_X1_dS_WyIEyFE
R2adj_X1_dSM_WyIEyFE
R2adj_X1_dS_WyIEn
R2adj_X1_dSM_WyIEn
R2adj_X1_dS_WyIEnFE
R2adj_X1_dSM_WyIEnFE

#Residual SE:
sigma(NLS_X1_dS_WyIEy)
sigma(NLS_X1_dSM_WyIEy)
sigma(NLS_X1_dS_WyIEyFE)
sigma(NLS_X1_dSM_WyIEyFE)
sigma(NLS_X1_dS_WyIEn)
sigma(NLS_X1_dSM_WyIEn)
sigma(NLS_X1_dS_WyIEnFE)
sigma(NLS_X1_dSM_WyIEnFE)


#Table A.11: RC min parties FTCU:----

#create data for Minority parties:
dM_WyIEy=dG_#base data from dG_

#create merging variable with unique codingL
dM_WyIEy$mer=paste0(dM_WyIEy$RILEg,dM_WyIEy$RILEp1,dM_WyIEy$RILEp2,dM_WyIEy$RILEp3,dM_WyIEy$RILEp4,dM_WyIEy$RILEp5
                    ,dM_WyIEy$RILEp6,dM_WyIEy$RILEp7,dM_WyIEy$RILEp8,dM_WyIEy$RILEp9,dM_WyIEy$RILEp10)
which(duplicated(dM_WyIEy$mer))

dt=dU_WyIEy[,c("RILEg",paste0("RILEp", 1:20),paste0("GOVp", 1:20,"_S"))]#base data for supporting parties
dt$mer=paste0(dt$RILEg,dt$RILEp1,dt$RILEp2,dt$RILEp3,dt$RILEp4,dt$RILEp5
                    ,dt$RILEp6,dt$RILEp7,dt$RILEp8,dt$RILEp9,dt$RILEp10)
which(duplicated(dt$mer))
dt=dt[,!(names(dt) %in% c("RILEg",paste0("RILEp", 1:20)))]

dM_WyIEy=merge(dM_WyIEy,dt,by=c("mer"),all.x=T,all.y=F)#merge
dM_WyIEy[,c(paste0("GOVp", 1:20,"_S"))][is.na(dM_WyIEy[,c(paste0("GOVp", 1:20,"_S"))])]=0#Na to 0 for data analysis

dM_WyIEy$GOVm1=ifelse(((dM_WyIEy$GOVp1==0)&(dM_WyIEy$GOVp1_S==0)&!is.na(dM_WyIEy$RILEp1)),1,0)#dummy being in the minority
dM_WyIEy$GOVm2=ifelse(((dM_WyIEy$GOVp2==0)&(dM_WyIEy$GOVp2_S==0)&!is.na(dM_WyIEy$RILEp2)),1,0)
dM_WyIEy$GOVm3=ifelse(((dM_WyIEy$GOVp3==0)&(dM_WyIEy$GOVp3_S==0)&!is.na(dM_WyIEy$RILEp3)),1,0)
dM_WyIEy$GOVm4=ifelse(((dM_WyIEy$GOVp4==0)&(dM_WyIEy$GOVp4_S==0)&!is.na(dM_WyIEy$RILEp4)),1,0)
dM_WyIEy$GOVm5=ifelse(((dM_WyIEy$GOVp5==0)&(dM_WyIEy$GOVp5_S==0)&!is.na(dM_WyIEy$RILEp5)),1,0)
dM_WyIEy$GOVm6=ifelse(((dM_WyIEy$GOVp6==0)&(dM_WyIEy$GOVp6_S==0)&!is.na(dM_WyIEy$RILEp6)),1,0)
dM_WyIEy$GOVm7=ifelse(((dM_WyIEy$GOVp7==0)&(dM_WyIEy$GOVp7_S==0)&!is.na(dM_WyIEy$RILEp7)),1,0)
dM_WyIEy$GOVm8=ifelse(((dM_WyIEy$GOVp8==0)&(dM_WyIEy$GOVp8_S==0)&!is.na(dM_WyIEy$RILEp8)),1,0)
dM_WyIEy$GOVm9=ifelse(((dM_WyIEy$GOVp9==0)&(dM_WyIEy$GOVp9_S==0)&!is.na(dM_WyIEy$RILEp9)),1,0)
dM_WyIEy$GOVm10=ifelse(((dM_WyIEy$GOVp10==0)&(dM_WyIEy$GOVp10_S==0)&!is.na(dM_WyIEy$RILEp10)),1,0)
dM_WyIEy$GOVm11=ifelse(((dM_WyIEy$GOVp11==0)&(dM_WyIEy$GOVp11_S==0)&!is.na(dM_WyIEy$RILEp11)),1,0)
dM_WyIEy$GOVm12=ifelse(((dM_WyIEy$GOVp12==0)&(dM_WyIEy$GOVp12_S==0)&!is.na(dM_WyIEy$RILEp12)),1,0)
dM_WyIEy$GOVm13=ifelse(((dM_WyIEy$GOVp13==0)&(dM_WyIEy$GOVp13_S==0)&!is.na(dM_WyIEy$RILEp13)),1,0)
dM_WyIEy$GOVm14=ifelse(((dM_WyIEy$GOVp14==0)&(dM_WyIEy$GOVp14_S==0)&!is.na(dM_WyIEy$RILEp14)),1,0)
dM_WyIEy$GOVm15=ifelse(((dM_WyIEy$GOVp15==0)&(dM_WyIEy$GOVp15_S==0)&!is.na(dM_WyIEy$RILEp15)),1,0)
dM_WyIEy$GOVm16=ifelse(((dM_WyIEy$GOVp16==0)&(dM_WyIEy$GOVp16_S==0)&!is.na(dM_WyIEy$RILEp16)),1,0)
dM_WyIEy$GOVm17=ifelse(((dM_WyIEy$GOVp17==0)&(dM_WyIEy$GOVp17_S==0)&!is.na(dM_WyIEy$RILEp17)),1,0)
dM_WyIEy$GOVm18=ifelse(((dM_WyIEy$GOVp18==0)&(dM_WyIEy$GOVp18_S==0)&!is.na(dM_WyIEy$RILEp18)),1,0)
dM_WyIEy$GOVm19=ifelse(((dM_WyIEy$GOVp19==0)&(dM_WyIEy$GOVp19_S==0)&!is.na(dM_WyIEy$RILEp19)),1,0)
dM_WyIEy$GOVm20=ifelse(((dM_WyIEy$GOVp20==0)&(dM_WyIEy$GOVp20_S==0)&!is.na(dM_WyIEy$RILEp20)),1,0)

dM_WyIEy$st=rowSums(dM_WyIEy[, names(dM_WyIEy) %in% c(paste0("SEATSp",1:20))])#tot seats (coded)
dM_WyIEy$SEATShPLp1=log(dM_WyIEy$SEATSp1/dM_WyIEy$st)#log seat share in parliament (coded)
dM_WyIEy$SEATShPLp2=log(dM_WyIEy$SEATSp2/dM_WyIEy$st)
dM_WyIEy$SEATShPLp3=log(dM_WyIEy$SEATSp3/dM_WyIEy$st)
dM_WyIEy$SEATShPLp4=log(dM_WyIEy$SEATSp4/dM_WyIEy$st)
dM_WyIEy$SEATShPLp5=log(dM_WyIEy$SEATSp5/dM_WyIEy$st)
dM_WyIEy$SEATShPLp6=log(dM_WyIEy$SEATSp6/dM_WyIEy$st)
dM_WyIEy$SEATShPLp7=log(dM_WyIEy$SEATSp7/dM_WyIEy$st)
dM_WyIEy$SEATShPLp8=log(dM_WyIEy$SEATSp8/dM_WyIEy$st)
dM_WyIEy$SEATShPLp9=log(dM_WyIEy$SEATSp9/dM_WyIEy$st)
dM_WyIEy$SEATShPLp10=log(dM_WyIEy$SEATSp10/dM_WyIEy$st)
dM_WyIEy$SEATShPLp11=log(dM_WyIEy$SEATSp11/dM_WyIEy$st)
dM_WyIEy$SEATShPLp12=log(dM_WyIEy$SEATSp12/dM_WyIEy$st)
dM_WyIEy$SEATShPLp13=log(dM_WyIEy$SEATSp13/dM_WyIEy$st)
dM_WyIEy$SEATShPLp14=log(dM_WyIEy$SEATSp14/dM_WyIEy$st)
dM_WyIEy$SEATShPLp15=log(dM_WyIEy$SEATSp15/dM_WyIEy$st)
dM_WyIEy$SEATShPLp16=log(dM_WyIEy$SEATSp16/dM_WyIEy$st)
dM_WyIEy$SEATShPLp17=log(dM_WyIEy$SEATSp17/dM_WyIEy$st)
dM_WyIEy$SEATShPLp18=log(dM_WyIEy$SEATSp18/dM_WyIEy$st)
dM_WyIEy$SEATShPLp19=log(dM_WyIEy$SEATSp19/dM_WyIEy$st)
dM_WyIEy$SEATShPLp20=log(dM_WyIEy$SEATSp20/dM_WyIEy$st)

dM_WyIEn=dM_WyIEy[dM_WyIEy$AfterElection==1,]#data with only Governments After Elections

#create matrices to pass to function:
X1columns_g=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVp",1:20),paste0("SEATShPLp",1:20))
X1columns_s=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVp",1:20,"_S"),paste0("SEATShPLp",1:20))
X1columns_m=c(paste0("RILEg"),paste0("RILEp",1:20),paste0("GOVm",1:20),paste0("SEATShPLp",1:20))

X1_dM_WyIEy_g=dM_WyIEy[,c(X1columns_g)]
X1_dM_WyIEy_g[is.na(X1_dM_WyIEy_g)]=-99
X1_dM_WyIEy_g=as.matrix(X1_dM_WyIEy_g)
X1_dM_WyIEy_s=dM_WyIEy[,c(X1columns_s)]
X1_dM_WyIEy_s[is.na(X1_dM_WyIEy_s)]=-99
X1_dM_WyIEy_s=as.matrix(X1_dM_WyIEy_s)
X1_dM_WyIEy_m=dM_WyIEy[,c(X1columns_m)]
X1_dM_WyIEy_m[is.na(X1_dM_WyIEy_m)]=-99
X1_dM_WyIEy_m=as.matrix(X1_dM_WyIEy_m)

X1_dM_WyIEn_g=dM_WyIEn[,c(X1columns_g)]
X1_dM_WyIEn_g[is.na(X1_dM_WyIEn_g)]=-99
X1_dM_WyIEn_g=as.matrix(X1_dM_WyIEn_g)
X1_dM_WyIEn_s=dM_WyIEn[,c(X1columns_s)]
X1_dM_WyIEn_s[is.na(X1_dM_WyIEn_s)]=-99
X1_dM_WyIEn_s=as.matrix(X1_dM_WyIEn_s)
X1_dM_WyIEn_m=dM_WyIEn[,c(X1columns_m)]
X1_dM_WyIEn_m[is.na(X1_dM_WyIEn_m)]=-99
X1_dM_WyIEn_m=as.matrix(X1_dM_WyIEn_m)

FE_dM_WyIEy=FE_dG(dM_WyIEy)#creating FE matrices
FE_dM_WyIEn=FE_dG(dM_WyIEn)

#regressions:
NLS_X1_dS_WyIEy=nls(RILEg~NLSc_dFSU1ii(b0, bgs, bng, Xg=X1_dM_WyIEy_g, Xs=X1_dM_WyIEy_s)
                    ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005))
                    ,data=dM_WyIEy)#
summary(NLS_X1_dS_WyIEy)
nrow(dM_WyIEy)
R2_X1_dS_WyIEy=1-sum((residuals(NLS_X1_dS_WyIEy))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dS_WyIEy
R2adj_X1_dS_WyIEy=1-(1-R2_X1_dS_WyIEy)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEy)))
R2adj_X1_dS_WyIEy
NLS_X1_dSM_WyIEy=nls(RILEg~NLSc_dFSU1iii(b0,bgs,bng,bm,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,Xm=X1_dM_WyIEy_m)
                     ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005))
                     ,data=dM_WyIEy)#
summary(NLS_X1_dSM_WyIEy)
R2_X1_dSM_WyIEy=1-sum((residuals(NLS_X1_dSM_WyIEy))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dSM_WyIEy
R2adj_X1_dSM_WyIEy=1-(1-R2_X1_dSM_WyIEy)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEy)))
R2adj_X1_dSM_WyIEy

NLS_X1_dS_WyIEyFE=nls(RILEg~NLSfe_dFSU1ii(bfe,bgs,bng,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,FE=FE_dM_WyIEy)
                      ,start = list(bgs=c(0.0005),bng=c(0.0005)
                                    ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                      ,data=dM_WyIEy)
summary(NLS_X1_dS_WyIEyFE)
R2_X1_dS_WyIEyFE=1-sum((residuals(NLS_X1_dS_WyIEyFE))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dS_WyIEyFE
R2adj_X1_dS_WyIEyFE=1-(1-R2_X1_dS_WyIEyFE)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEyFE)))
R2adj_X1_dS_WyIEyFE
NLS_X1_dSM_WyIEyFE=nls(RILEg~NLSfe_dFSU1iii(bfe,bgs,bng,bm,Xg=X1_dM_WyIEy_g,Xs=X1_dM_WyIEy_s,Xm=X1_dM_WyIEy_m,FE=FE_dM_WyIEy)
                       ,start = list(bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005)
                                     ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                       ,data=dM_WyIEy)
summary(NLS_X1_dSM_WyIEyFE)
R2_X1_dSM_WyIEyFE=1-sum((residuals(NLS_X1_dSM_WyIEyFE))^2)/sum((dM_WyIEy$RILEg-mean(dM_WyIEy$RILEg))^2)
R2_X1_dSM_WyIEyFE
R2adj_X1_dSM_WyIEyFE=1-(1-R2_X1_dSM_WyIEyFE)*((length(dM_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEyFE)))
R2adj_X1_dSM_WyIEyFE

NLS_X1_dS_WyIEn=nls(RILEg~NLSc_dFSU1ii(b0, bgs, bng, Xg=X1_dM_WyIEn_g, Xs=X1_dM_WyIEn_s)
                    ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005))
                    ,data=dM_WyIEn)#
summary(NLS_X1_dS_WyIEn)
R2_X1_dS_WyIEn=1-sum((residuals(NLS_X1_dS_WyIEn))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dS_WyIEn
R2adj_X1_dS_WyIEn=1-(1-R2_X1_dS_WyIEn)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEn)))
R2adj_X1_dS_WyIEn
NLS_X1_dSM_WyIEn=nls(RILEg~NLSc_dFSU1iii(b0,bgs,bng,bm,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,Xm=X1_dM_WyIEn_m)
                     ,start = list(b0=1,bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005))
                     ,data=dM_WyIEn)#
summary(NLS_X1_dSM_WyIEn)
R2_X1_dSM_WyIEn=1-sum((residuals(NLS_X1_dSM_WyIEn))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dSM_WyIEn
R2adj_X1_dSM_WyIEn=1-(1-R2_X1_dSM_WyIEn)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEn)))
R2adj_X1_dSM_WyIEn

NLS_X1_dS_WyIEnFE=nls(RILEg~NLSfe_dFSU1ii(bfe,bgs,bng,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,FE=FE_dM_WyIEn)
                      ,start = list(bgs=c(0.0005),bng=c(0.0005)
                                    ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                      ,data=dM_WyIEn)
summary(NLS_X1_dS_WyIEnFE)
R2_X1_dS_WyIEnFE=1-sum((residuals(NLS_X1_dS_WyIEnFE))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dS_WyIEnFE
R2adj_X1_dS_WyIEnFE=1-(1-R2_X1_dS_WyIEnFE)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dS_WyIEnFE)))
R2adj_X1_dS_WyIEnFE
NLS_X1_dSM_WyIEnFE=nls(RILEg~NLSfe_dFSU1iii(bfe,bgs,bng,bm,Xg=X1_dM_WyIEn_g,Xs=X1_dM_WyIEn_s,Xm=X1_dM_WyIEn_m,FE=FE_dM_WyIEn)
                       ,start = list(bgs=c(0.0005),bng=c(0.0005),bm=c(0.0005)
                                     ,bfe=c(0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005,0.0005))
                       ,data=dM_WyIEn)
summary(NLS_X1_dSM_WyIEnFE)
R2_X1_dSM_WyIEnFE=1-sum((residuals(NLS_X1_dSM_WyIEnFE))^2)/sum((dM_WyIEn$RILEg-mean(dM_WyIEn$RILEg))^2)
R2_X1_dSM_WyIEnFE
R2adj_X1_dSM_WyIEnFE=1-(1-R2_X1_dSM_WyIEnFE)*((length(dM_WyIEn$RILEg)-1)/(df.residual(NLS_X1_dSM_WyIEnFE)))
R2adj_X1_dSM_WyIEnFE

#exporting regressions to latex:
x<- matrix(rnorm(nrow(dG_WyIEy)*13), ncol = 13)
x<- data.frame(x)
colnames(x)<- c("RILEg", "bgs", "bng", "bm"
                ,"FE1", "FE2", "FE3", "FE4", "FE5", "FE6", "FE7", "FE8", "FE9")
mod1<- lm(RILEg~ bgs+bng, data = x)
coefs1<- NLS_X1_dS_WyIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_dS_WyIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ bgs+bng+bm, data = x)
coefs2<- NLS_X1_dSM_WyIEy$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X1_dSM_WyIEy)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ bgs+bng+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs3<- NLS_X1_dS_WyIEyFE$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X1_dS_WyIEyFE)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ bgs+bng+bm+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs4<- NLS_X1_dSM_WyIEyFE$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X1_dSM_WyIEyFE)))
names(se4)<- names(coefs4)
mod5<- lm(RILEg~ bgs+bng, data = x)
coefs5<- NLS_X1_dS_WyIEn$m$getPars()
names(coefs5)<- names(mod5$coefficients)
se5<- sqrt(diag(vcov(NLS_X1_dS_WyIEn)))
names(se5)<- names(coefs5)
mod6<- lm(RILEg~ bgs+bng+bm, data = x)
coefs6<- NLS_X1_dSM_WyIEn$m$getPars()
names(coefs6)<- names(mod6$coefficients)
se6<- sqrt(diag(vcov(NLS_X1_dSM_WyIEn)))
names(se6)<- names(coefs6)
mod7<- lm(RILEg~ bgs+bng+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs7<- NLS_X1_dS_WyIEnFE$m$getPars()
names(coefs7)<- names(mod7$coefficients)
se7<- sqrt(diag(vcov(NLS_X1_dS_WyIEnFE)))
names(se7)<- names(coefs7)
mod8<- lm(RILEg~ bgs+bng+bm+ FE1 + FE2 + FE3 + FE4 + FE5 + FE6 + FE7 + FE8 -1, data = x)
coefs8<- NLS_X1_dSM_WyIEnFE$m$getPars()
names(coefs8)<- names(mod8$coefficients)
se8<- sqrt(diag(vcov(NLS_X1_dSM_WyIEnFE)))
names(se8)<- names(coefs8)
stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,
          coef = list(coefs1,coefs2,coefs3,coefs4,coefs5,coefs6,coefs7,coefs8),
          se = list(se1,se2,se3,se4,se5,se6,se7,se8),
          covariate.labels = c("Log Seat Share", "Supporting Parties", "Opposition Parties"),
          type = "latex",
          out = "TableA.11.tex")

#Pseudo R2:
R2_X1_dS_WyIEy
R2_X1_dSM_WyIEy
R2_X1_dS_WyIEyFE
R2_X1_dSM_WyIEyFE
R2_X1_dS_WyIEn
R2_X1_dSM_WyIEn
R2_X1_dS_WyIEnFE
R2_X1_dSM_WyIEnFE

#Pseudo Adjusted R2:
R2adj_X1_dS_WyIEy
R2adj_X1_dSM_WyIEy
R2adj_X1_dS_WyIEyFE
R2adj_X1_dSM_WyIEyFE
R2adj_X1_dS_WyIEn
R2adj_X1_dSM_WyIEn
R2adj_X1_dS_WyIEnFE
R2adj_X1_dSM_WyIEnFE

#Residual SE:
sigma(NLS_X1_dS_WyIEy)
sigma(NLS_X1_dSM_WyIEy)
sigma(NLS_X1_dS_WyIEyFE)
sigma(NLS_X1_dSM_WyIEyFE)
sigma(NLS_X1_dS_WyIEn)
sigma(NLS_X1_dSM_WyIEn)
sigma(NLS_X1_dS_WyIEnFE)
sigma(NLS_X1_dSM_WyIEnFE)


#Tables A.4-A.6: RC power indexes and MIW:----

##Table A.6 model 1 (run again main model from Table 2 as comparison model):
X3_134columns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)#create data:
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
X3_134=dG_WyIEy[, names(dG_WyIEy) %in% X3_134columns]
X3_134[is.na(X3_134)]=-99
X3_134=as.matrix(X3_134)
#regressions:
NLS_X3_134c=nls(RILEg~NLSc(b0, b, X=X3_134)
                ,start = list(b0=1,b=c(0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X3_134c)


##Table A.4 model 1 (ISmodel discrimination):
XIScolumns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("ISp", 1:20))
XIS_dG_WyIEy=dG_WyIEy[, names(dG_WyIEy) %in% XIScolumns]
XIS_dG_WyIEy[is.na(XIS_dG_WyIEy)]=-99
XIS_dG_WyIEy=as.matrix(XIS_dG_WyIEy)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_XIS_dG_WyIEy=nls(RILEg~NLSc(b0, b, X=XIS_dG_WyIEy)
                     ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEy)
summary(NLS_XIS_dG_WyIEy)
R2_XIS_dG_WyIEy=1-sum((residuals(NLS_XIS_dG_WyIEy))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_XIS_dG_WyIEy
R2adj_XIS_dG_WyIEy=1-(1-R2_XIS_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_XIS_dG_WyIEy)))
R2adj_XIS_dG_WyIEy
#test:
vuongtest(NLS_XIS_dG_WyIEy,NLS_X3_134c)
icci(NLS_XIS_dG_WyIEy,NLS_X3_134c)


##Table A.5 model 1 (same model in Table 1 model 1):
summary(NLS_X1_dG_WyIEy)
R2_X1_dG_WyIEy=1-sum((residuals(NLS_X1_dG_WyIEy))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1_dG_WyIEy
R2adj_X1_dG_WyIEy=1-(1-R2_X1_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X1_dG_WyIEy)))
R2adj_X1_dG_WyIEy
#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_dG_WyIEy)$coefficients[2,4]
2*pt(abs((NLS_X1_dG_WyIEy$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEy)))[2]),summary(NLS_X1_dG_WyIEy)$df[2],lower.tail = FALSE)


##Table A.5 model 2 (IS same model):
X2_1IScolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ISp", 1:20))
X2_1IS=dG_WyIEy[, names(dG_WyIEy) %in% X2_1IScolumns]
X2_1IS[is.na(X2_1IS)]=-99
X2_1IS=as.matrix(X2_1IS)
X4_134IScolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                  , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20), paste0("ISp", 1:20))
X4_134IS=dG_WyIEy[, names(dG_WyIEy) %in% X4_134IScolumns]
X4_134IS[is.na(X4_134IS)]=-99
X4_134IS=as.matrix(X4_134IS)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_X2_1ISc=nls(RILEg~NLSc(b0, b, X=X2_1IS)
                ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_1ISc)
R2_X1ISc_dG_WyIEy=1-sum((residuals(NLS_X2_1ISc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1ISc_dG_WyIEy
R2adj_X1ISc_dG_WnIEy=1-(1-R2_X1ISc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_1ISc)))
R2adj_X1ISc_dG_WnIEy
summary(NLS_X2_1ISc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X2_1ISc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1ISc)))[2])
     ,summary(NLS_X2_1ISc)$df[2],lower.tail = FALSE)#p-value null beta=1


##Table A.6 model 2 (IS same model):
NLS_X4_134ISc=nls(RILEg~NLSc(b0, b, X=X4_134IS)
                  ,start = list(b0=1,b=c(0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X4_134ISc)
R2_X4134ISc_dG_WyIEy=1-sum((residuals(NLS_X4_134ISc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X4134ISc_dG_WyIEy
R2adj_X4134ISc_dG_WnIEy=1-(1-R2_X4134ISc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X4_134ISc)))
R2adj_X4134ISc_dG_WnIEy
summary(NLS_X4_134ISc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X4_134ISc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134ISc)))[2])
     ,summary(NLS_X4_134ISc)$df[2],lower.tail = FALSE)#p-value null beta=1


##Table A.4 model 2 (IBmodel discrimination):
XIBcolumns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("IBp", 1:20))
XIB_dG_WyIEy=dG_WyIEy[, names(dG_WyIEy) %in% XIBcolumns]
XIB_dG_WyIEy[is.na(XIB_dG_WyIEy)]=-99
XIB_dG_WyIEy=as.matrix(XIB_dG_WyIEy)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_XIB_dG_WyIEy=nls(RILEg~NLSc(b0, b, X=XIB_dG_WyIEy)
                     ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEy)
summary(NLS_XIB_dG_WyIEy)
R2_XIB_dG_WyIEy=1-sum((residuals(NLS_XIB_dG_WyIEy))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_XIB_dG_WyIEy
R2adj_XIB_dG_WyIEy=1-(1-R2_XIB_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_XIB_dG_WyIEy)))
R2adj_XIB_dG_WyIEy
#test:
vuongtest(NLS_XIB_dG_WyIEy,NLS_X3_134c)
icci(NLS_XIB_dG_WyIEy,NLS_X3_134c)


##Table A.5 model 3 (IB same model):
X2_1IBcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("IBp", 1:20))
X2_1IB=dG_WyIEy[, names(dG_WyIEy) %in% X2_1IBcolumns]
X2_1IB[is.na(X2_1IB)]=-99
X2_1IB=as.matrix(X2_1IB)
X4_134IBcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                  , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20), paste0("IBp", 1:20))
X4_134IB=dG_WyIEy[, names(dG_WyIEy) %in% X4_134IBcolumns]
X4_134IB[is.na(X4_134IB)]=-99
X4_134IB=as.matrix(X4_134IB)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_X2_1IBc=nls(RILEg~NLSc(b0, b, X=X2_1IB)
                ,start = list(b0=1,b=c(0.0005,0.0005))
                ,nls.control(maxiter = 1000, minFactor = 0.00000001, warnOnly = T)
                ,data=dG_WyIEy)
summary(NLS_X2_1IBc)
R2_X1IBc_dG_WyIEy=1-sum((residuals(NLS_X2_1IBc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1IBc_dG_WyIEy
R2adj_X1IBc_dG_WnIEy=1-(1-R2_X1IBc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_1IBc)))
R2adj_X1IBc_dG_WnIEy
summary(NLS_X2_1IBc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X2_1IBc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1IBc)))[2])
     ,summary(NLS_X2_1IBc)$df[2],lower.tail = FALSE)#p-value null beta=1


##Table A.6 model 3 (IB same model):
NLS_X4_134IBc=nls(RILEg~NLSc(b0, b, X=X4_134IB)
                  ,start = list(b0=1,b=c(0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X4_134IBc)
R2_X4134IBc_dG_WyIEy=1-sum((residuals(NLS_X4_134IBc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X4134IBc_dG_WyIEy
R2adj_X4134IBc_dG_WnIEy=1-(1-R2_X4134IBc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X4_134IBc)))
R2adj_X4134IBc_dG_WnIEy
summary(NLS_X4_134IBc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X4_134IBc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134IBc)))[2])
     ,summary(NLS_X4_134IBc)$df[2],lower.tail = FALSE)#p-value null beta=1


##Table A.4 model 3 (MIWmodel discrimination):
XMIWcolumns=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20), paste0("MIWp", 1:20))
XMIW_dG_WyIEy=dG_WyIEy[, names(dG_WyIEy) %in% XMIWcolumns]
XMIW_dG_WyIEy[is.na(XMIW_dG_WyIEy)]=-99
XMIW_dG_WyIEy=as.matrix(XMIW_dG_WyIEy)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_XMIW_dG_WyIEy=nls(RILEg~NLSc(b0, b, X=XMIW_dG_WyIEy)
                      ,start = list(b0=1,b=c(0.0005)),data=dG_WyIEy)
summary(NLS_XMIW_dG_WyIEy)
R2_XMIW_dG_WyIEy=1-sum((residuals(NLS_XMIW_dG_WyIEy))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_XMIW_dG_WyIEy
R2adj_XMIW_dG_WyIEy=1-(1-R2_XMIW_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_XMIW_dG_WyIEy)))
R2adj_XMIW_dG_WyIEy
#test:
vuongtest(NLS_XMIW_dG_WyIEy,NLS_X3_134c)
icci(NLS_XMIW_dG_WyIEy,NLS_X3_134c)


##Table A.5 model 4 (MIW same model):
X2_1MIWcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                 , paste0("SEATShCLp", 1:20), paste0("MIWp", 1:20))
X2_1MIW=dG_WyIEy[, names(dG_WyIEy) %in% X2_1MIWcolumns]
X2_1MIW[is.na(X2_1MIW)]=-99
X2_1MIW=as.matrix(X2_1MIW)
X4_134MIWcolumns=c(paste0("RILEg"), paste0("GOVp", 1:20), paste0("RILEp", 1:20)
                   , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20), paste0("MIWp", 1:20))
X4_134MIW=dG_WyIEy[, names(dG_WyIEy) %in% X4_134MIWcolumns]
X4_134MIW[is.na(X4_134MIW)]=-99
X4_134MIW=as.matrix(X4_134MIW)
FE_dG_WyIEy=FE_dG(dG_WyIEy)#creating FE matrices
#regressions:
NLS_X2_1MIWc=nls(RILEg~NLSc(b0, b, X=X2_1MIW)
                 ,start = list(b0=1,b=c(0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X2_1MIWc)
R2_X1MIWc_dG_WyIEy=1-sum((residuals(NLS_X2_1MIWc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X1MIWc_dG_WyIEy
R2adj_X1MIWc_dG_WnIEy=1-(1-R2_X1MIWc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X2_1MIWc)))
R2adj_X1MIWc_dG_WnIEy
summary(NLS_X2_1MIWc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X2_1MIWc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1MIWc)))[2])
     ,summary(NLS_X2_1MIWc)$df[2],lower.tail = FALSE)#p-value null beta=1


##Table A.6 model 4 (MIW same model):
NLS_X4_134MIWc=nls(RILEg~NLSc(b0, b, X=X4_134MIW)
                   ,start = list(b0=1,b=c(0.0005,0.0005,0.0005,0.0005)),data=dG_WyIEy)
summary(NLS_X4_134MIWc)
R2_X4134MIWc_dG_WyIEy=1-sum((residuals(NLS_X4_134MIWc))^2)/sum((dG_WyIEy$RILEg-mean(dG_WyIEy$RILEg))^2)
R2_X4134MIWc_dG_WyIEy
R2adj_X4134MIWc_dG_WnIEy=1-(1-R2_X4134MIWc_dG_WyIEy)*((length(dG_WyIEy$RILEg)-1)/(df.residual(NLS_X4_134MIWc)))
R2adj_X4134MIWc_dG_WnIEy
summary(NLS_X4_134MIWc)$coefficients[2,4]#p-value null beta=0
2*pt(abs((NLS_X4_134MIWc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134MIWc)))[2])
     ,summary(NLS_X4_134MIWc)$df[2],lower.tail = FALSE)#p-value null beta=1


##exporting regressions to latex - TableA.4:
x<- matrix(rnorm(nrow(dG_WyIEy)*13), ncol = 13)
x<- data.frame(x)
colnames(x)<- c("RILEg", "bssi", "bbi", "bmiw")
mod1<- lm(RILEg~ bssi, data = x)
coefs1<- NLS_XIS_dG_WyIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_XIS_dG_WyIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ bbi, data = x)
coefs2<- NLS_XIB_dG_WyIEy$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_XIB_dG_WyIEy)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ bmiw, data = x)
coefs3<- NLS_XMIW_dG_WyIEy$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_XMIW_dG_WyIEy)))
names(se3)<- names(coefs3)
stargazer(mod1,mod2,mod3,
          coef = list(coefs1,coefs2,coefs3),
          se = list(se1,se2,se3),
          covariate.labels = c("Shapley-Shubik Index", "Banzhaf Index", "Minimum Integer Weights"),
          type = "latex",
          out = "TableA.4.tex")

#Pseudo R2:
R2_XIS_dG_WyIEy
R2_XIB_dG_WyIEy
R2_XMIW_dG_WyIEy

#Pseudo Adjusted R2:
R2adj_XIS_dG_WyIEy
R2adj_XIB_dG_WyIEy
R2adj_XMIW_dG_WyIEy

#Residual SE:
sigma(NLS_XIS_dG_WyIEy)
sigma(NLS_XIB_dG_WyIEy)
sigma(NLS_XMIW_dG_WyIEy)

#Vuong Test:
vuongtest(NLS_XIS_dG_WyIEy,NLS_X3_134c)$LRTstat
vuongtest(NLS_XIS_dG_WyIEy,NLS_X3_134c)$p_LRT$B
vuongtest(NLS_XIB_dG_WyIEy,NLS_X3_134c)$LRTstat
vuongtest(NLS_XIB_dG_WyIEy,NLS_X3_134c)$p_LRT$B
vuongtest(NLS_XMIW_dG_WyIEy,NLS_X3_134c)$LRTstat
vuongtest(NLS_XMIW_dG_WyIEy,NLS_X3_134c)$p_LRT$B

#Akaike Information Criterion:
icci(NLS_XIS_dG_WyIEy,NLS_X3_134c)$AIC[1]
icci(NLS_XIB_dG_WyIEy,NLS_X3_134c)$AIC[1]
icci(NLS_XMIW_dG_WyIEy,NLS_X3_134c)$AIC[1]

#Bayesian Information Criterion:
icci(NLS_XIS_dG_WyIEy,NLS_X3_134c)$BIC[1]
icci(NLS_XIB_dG_WyIEy,NLS_X3_134c)$BIC[1]
icci(NLS_XMIW_dG_WyIEy,NLS_X3_134c)$BIC[1]


##exporting regressions to latex - TableA.5:
x<- matrix(rnorm(nrow(dG_WyIEy)*5), ncol = 5)
x<- data.frame(x)
colnames(x)<- c("RILEg", "blss", "bssi", "bbi", "bmiw")
mod1<- lm(RILEg~ blss, data = x)
coefs1<- NLS_X1_dG_WyIEy$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X1_dG_WyIEy)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ blss+bssi, data = x)
coefs2<- NLS_X2_1ISc$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X2_1ISc)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ blss+bbi, data = x)
coefs3<- NLS_X2_1IBc$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X2_1IBc)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ blss+bmiw, data = x)
coefs4<- NLS_X2_1MIWc$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X2_1MIWc)))
names(se4)<- names(coefs4)
stargazer(mod1, mod2, mod3, mod4,
          coef = list(coefs1, coefs2, coefs3, coefs4),
          se = list(se1, se2, se3, se4),
          covariate.labels = c("Log Seat Share", "Shapley-Shubik Index", "Banzhaf Index", "Minimum Integer Weights"),
          type = "latex",
          out = "TableA.5.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X1_dG_WyIEy)$coefficients[2,4]
summary(NLS_X2_1ISc)$coefficients[2,4]
summary(NLS_X2_1IBc)$coefficients[2,4]
summary(NLS_X2_1MIWc)$coefficients[2,4]
2*pt(abs((NLS_X1_dG_WyIEy$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X1_dG_WyIEy)))[2]),summary(NLS_X1_dG_WyIEy)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_1ISc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1ISc)))[2])
     ,summary(NLS_X2_1ISc)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_1IBc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1IBc)))[2])
     ,summary(NLS_X2_1IBc)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X2_1MIWc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X2_1MIWc)))[2])
     ,summary(NLS_X2_1MIWc)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X1_dG_WyIEy
R2_X1ISc_dG_WyIEy
R2_X1IBc_dG_WyIEy
R2_X1MIWc_dG_WyIEy

#Pseudo Adjusted R2:
R2adj_X1_dG_WyIEy
R2adj_X1ISc_dG_WnIEy
R2adj_X1IBc_dG_WnIEy
R2adj_X1MIWc_dG_WnIEy

#Residual SE:
sigma(NLS_X1_dG_WyIEy)
sigma(NLS_X2_1ISc)
sigma(NLS_X2_1IBc)
sigma(NLS_X2_1MIWc)


##exporting regressions to latex - TableA.6:
x<- matrix(rnorm(nrow(dG_WyIEy)*7), ncol = 7)
x<- data.frame(x)
colnames(x)<- c("RILEg", "blss", "b4", "b5", "bssi", "bbi", "bmiw")
mod1<- lm(RILEg~ blss+b4+b5, data = x)
coefs1<- NLS_X3_134c$m$getPars()
names(coefs1)<- names(mod1$coefficients)
se1<- sqrt(diag(vcov(NLS_X3_134c)))
names(se1)<- names(coefs1)
mod2<- lm(RILEg~ blss+b4+b5+bssi, data = x)
coefs2<- NLS_X4_134ISc$m$getPars()
names(coefs2)<- names(mod2$coefficients)
se2<- sqrt(diag(vcov(NLS_X4_134ISc)))
names(se2)<- names(coefs2)
mod3<- lm(RILEg~ blss+b4+b5+bbi, data = x)
coefs3<- NLS_X4_134IBc$m$getPars()
names(coefs3)<- names(mod3$coefficients)
se3<- sqrt(diag(vcov(NLS_X4_134IBc)))
names(se3)<- names(coefs3)
mod4<- lm(RILEg~ blss+b4+b5+bmiw, data = x)
coefs4<- NLS_X4_134MIWc$m$getPars()
names(coefs4)<- names(mod4$coefficients)
se4<- sqrt(diag(vcov(NLS_X4_134MIWc)))
names(se4)<- names(coefs4)
stargazer(mod1, mod2, mod3, mod4,
          coef = list(coefs1, coefs2, coefs3, coefs4),
          se = list(se1, se2, se3, se4),
          covariate.labels = c("Log Seat Share", "Abs. Distance to Median Legislative Party", "Formateur", "Shapley-Shubik Index", "Banzhaf Index", "Minimum Integer Weights"),
          type = "latex",
          out = "TableA.6.tex")

#table p-values test with null beta=0 and beta=1:
summary(NLS_X3_134c)$coefficients[2,4]
summary(NLS_X4_134ISc)$coefficients[2,4]
summary(NLS_X4_134IBc)$coefficients[2,4]
summary(NLS_X4_134MIWc)$coefficients[2,4]
2*pt(abs((NLS_X3_134c$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X3_134c)))[2]),summary(NLS_X3_134c)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X4_134ISc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134ISc)))[2])
     ,summary(NLS_X4_134ISc)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X4_134IBc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134IBc)))[2])
     ,summary(NLS_X4_134IBc)$df[2],lower.tail = FALSE)
2*pt(abs((NLS_X4_134MIWc$m$getPars()[2]-1)/sqrt(diag(vcov(NLS_X4_134MIWc)))[2])
     ,summary(NLS_X4_134MIWc)$df[2],lower.tail = FALSE)

#Pseudo R2:
R2_X3_134c
R2_X4134ISc_dG_WyIEy
R2_X4134IBc_dG_WyIEy
R2_X4134MIWc_dG_WyIEy

#Pseudo Adjusted R2:
R2adj_X3_134c
R2adj_X4134ISc_dG_WnIEy
R2adj_X4134IBc_dG_WnIEy
R2adj_X4134MIWc_dG_WnIEy

#Residual SE:
sigma(NLS_X3_134c)
sigma(NLS_X4_134ISc)
sigma(NLS_X4_134IBc)
sigma(NLS_X4_134MIWc)


#MARGINAL EFFECTS:----

#summary main dataset for paper:
summary(dG_$RILEg)
sd(dG_$RILEg)#15 (15.94348)
Pmed=dG_[,names(dG_) %in% paste0("RILEp", 1:20)]
Pmed=as.matrix(Pmed)
min(Pmed,na.rm=TRUE)
max(Pmed,na.rm=TRUE)
median(Pmed,na.rm=TRUE)
sd(Pmed,na.rm=TRUE)#20 (20.43498) 

##average seat share governing and non-governing parties:
G=as.matrix(dM_WyIEy[,c(paste0("GOVp",1:20))])*as.matrix(dM_WyIEy[,c(paste0("SEATShPLp",1:20))])
G[G==0|G==-Inf]=NA
G=exp(G)
mean(G,na.rm=T)#avg seat share gov parties: 0.2061329
G=as.matrix(1-dM_WyIEy[,c(paste0("GOVp",1:20))])*as.matrix(dM_WyIEy[,c(paste0("SEATShPLp",1:20))])
G[G==0|G==-Inf]=NA
G=exp(G)
mean(G,na.rm=T)#avg seat share non gov parties: 0.1299095
rm(G)

summary(NLS_X3_134c)$coefficients[1,1]#constant
summary(NLS_X3_134c)$coefficients[2,1]#log seat share
summary(NLS_X3_134c)$coefficients[3,1]#distance to the median
summary(NLS_X3_134c)$coefficients[4,1]#formateur

NLScME=function(b0, b, X){ #function with constant == NLSc w/ print(predict)
  K=((ncol(X)-41)/20)
  XX=list()
  SbX=matrix(0,nrow = nrow(X), ncol = ncol(X))
  for (k in 1:K) {
    XX[[k]]=X[1:nrow(X),(42+(k-1)*20):(41+k*20)]
    XX[[k]]=b[k]*XX[[k]]
  }
  SbX=Reduce('+', XX)
  eSbX=exp(SbX)
  eSbXG=eSbX*X[,22:41]
  eSbXGRsum=rowSums(eSbXG, na.rm = TRUE)
  eSbXGUnit=eSbXG/eSbXGRsum
  eSbXGUnitRILE=eSbXGUnit*X[,2:21]
  eSbXGUnitRILEsum=rowSums(eSbXGUnitRILE, na.rm = TRUE)
  predict=b0+eSbXGUnitRILEsum
  print(eSbXGUnit[,c(1:sum(X[1,22:41]==1,na.rm=T))])#print party weights
  print(predict)#print govt position
}

XME=matrix(0,nrow=3,ncol=101)
# colnames(XME)=colnames(X3_134)
colnames(XME)=c(paste0("RILEg"), paste0("RILEp", 1:20), paste0("GOVp", 1:20)
                , paste0("SEATShCLp", 1:20), paste0("ADMVPCp", 1:20), paste0("FORMp", 1:20))
XME[c(1,2),"RILEp1"]=-1*20#assign RILE
XME[c(3),"RILEp1"]=-12*20
XME[,"RILEp2"]=0
XME[,c("GOVp1","GOVp2")]=1#assign 2 coalition parties
XME[,c("SEATShCLp1","SEATShCLp2")]=log(0.5)#with equal log seat share (splitting total seat in the coalition)
XME[c(1,2),"ADMVPCp1"]=(1*20)/200#assign ADMVPCp
XME[c(3),"ADMVPCp1"]=(12*20)/200
XME[1,"FORMp1"]=1#assign FORMp
XME[2,"FORMp2"]=1
XME[3,"FORMp1"]=1
NLScME(b0=summary(NLS_X3_134c)$coefficients[1,1]
     ,b=c(summary(NLS_X3_134c)$coefficients[2,1]
          ,summary(NLS_X3_134c)$coefficients[3,1]
          ,summary(NLS_X3_134c)$coefficients[4,1])
     ,X=XME)
# [1,]  0.7855206  0.2144794
# [2,]  0.1320565  0.8679435
# [3,]  0.1280774  0.8719226
# [1]  -7.712701   5.356581 -22.740869
# -7.712701-5.356581=-13.06928/20~=2/3 govt position change switching formateur
# 0.86/0.21~=4 increase in weight for party 2 by becoming the formateur
# 1 to 12 sd move away party 1 to have same increase in weight for party 2



